arXiv: 1508.07903v2 [gr-qc] 13 Nov 2015 


Published in JCAP 11 (2015) 022 


The theory of stochastic cosmological 
lensing 


Pierre Fleury,“’^ Julien Larena,*^ Jean-Philippe Uzan“’^ 

“Institut d’Astrophysique de Paris, UMR 7095 du CNRS, 98 bis Bd Arago, 75014 Paris, 
France. 

^Sorbonne Universites, Institut Lagrange de Paris, 98 bis, Bd Arago, 75014 Paris, France. 
'^Department of Mathematics, Rhodes University, Grahamstown 6140, South Africa 

E-mail: fleury@iap.fr, j.larena@ru.ac.za, uzan@iap.fr 

Abstract. On the scale of the light beams subtended by small sources, e.g. supernovae, matter 
cannot be accurately described as a fluid, which questions the applicability of standard cosmic 
lensing to those cases. In this article, we propose a new formalism to deal with small-scale 
lensing as a diffusion process: the Sachs and Jacobi equations governing the propagation 
of narrow light beams are treated as Langevin equations. We derive the associated Fokker- 
Planck-Kolmogorov equations, and use them to deduce general analytical results on the mean 
and dispersion of the angular distance. This formalism is applied to random Einstein-Straus 
Swiss-cheese models, allowing us to: (1) show an explicit example of the involved calculations; 
(2) check the validity of the method against both ray-tracing simulations and direct numerical 
integration of the Langevin equation. As a byproduct, we obtain a post-Kantowski-Dyer- 
Roeder approximation, accounting for the effect of tidal distortions on the angular distance, 
in excellent agreement with numerical results. Besides, the dispersion of the angular distance 
is correctly reproduced in some regimes. 
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1 Introduction 

The understanding of light propagation in the Universe, in particular through the relation 
between distances and redshifts, is central for the interpretation of almost all cosmological 
observations. The standard approach consists in assuming that light propagates through a 
strictly homogeneous and isotropic Friedmann-Lemaitre (EL) spacetime [1], assumed to be a 
good model on cosmological scales.^ Such a crude—but surprisingly efficient—approximation 
can be refined by taking into account: (i) the actual non-comobility of both the light sources 
and the observer; (ii) the gravitational leasing caused by the large-scale structure. This more 
realistic description generally relies on the cosmological perturbation theory [5—7]. At first 
order, it essentially introduces a dispersion of the distance-redshift relation with respect to the 
background FL prediction [8-12], which can be partially corrected if a leasing map is known. 
There was recently an interesting debate on the bias potentially introdnced by second-order 
corrections: based on the calculations of Refs. [13, 14] (see also Refs. [15-17] for earlier 
results). Ref. [18] suggested that second-order leasing could significantly affect the standard 
interpretation of the cosmic microwave backgronnd (CMB) observations. Nevertheless, this 
statement tnrned ont to be inaccnrate, dne to confusions between several averaging schemes 
for the observable qnantities at stake [19-22]. 

This problem of determining the effect of inhomogeneities on light propagation can 
also been tackled in a nonperturbative way, e.g. by relying on toy models. The most 
common examples are Swiss-cheese models [23, 24], where inhomogeneities are introduced 
within a background FL spacetime by inserting spherical patches of another exact solntion of 
Einstein’s eqnation. Recent analyses generally exploit the Lemaitre-Tolman-Bondi (LTB) [25- 
39] or Szekeres [40-43] geometries as interior solutions, which aim at describing large-scale 
inhomogeneities snch as superclusters or cosmic voids (see also Refs. [44, 45]). Observations 
have also been connected to the cosmic coarse-graining and backreaction issues in the series 
of works [46-54] . 

^See however Refs. [2-4] for a recent debate on this specific issue. 
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All the above-mentioned approaches have in common that they describe matter in the 
Universe as a fluid. However, when it comes to narrow beams, such as those involved in 
supernova (SN) observations, this approximation should no longer hold.^ The applicability of 
the perturbation theory in this regime, in particular, has been questioned in Ref. [55]. This 
specific issue of how the dumpiness of the Universe affects the interpretation of cosmological 
observables was first raised by Zel’dovich [56] and Feynman [57]. The basic underlying idea is 
that in a clumpy medium, light mostly propagates through vacuum, and therefore experiences 
an underdense Universe. This stimulated a corpus of seminal articles [58-64], including the 
first analyses based on a Swiss-cheese model with Schwarzschild vacuoles [65-70]. Contrary 
to LTB or Szekeres holes, the latter aim at modelling relatively small gravitationally bound 
structures, such as individual galaxies or stars. The analysis of light propagation in such 
models resulted in the so-called Dyer-Roeder approximation—that we shall rather call the 
Kantowski-Dyer-Roeder (KDR) approximation in this article, the name of Kantowski being 
unfairly omitted in the literature. Its correspondence with Swiss-cheese models has been 
carefully rederived and numerically checked in Ref. [71], although its mathematical consistency 
was questioned in Refs. [49, 55]. Analyses based on other models than Swiss cheeses, albeit 
physically similar in the sense that they also describe universes made of point masses, have 
been proposed in Refs. [72-77]. When applied to the interpretation of SN data, these various 
approaches generically do And a bias in the measurement of the cosmological parameters, 
on the order of a few to more than ten percent [78-82]. It has been shown in Ref. [81] that 
such an effect improves the agreement between SN and CMB observations regarding the 
measurement of Hmo- 

While the KDR approximation may capture the main effects of the Universe’s dumpiness 
on the average distance-redshift relation, it does not tell anything about its dispersion, and 
a fortiori about its higher-order moments. Model-based approaches do not in principle 
suffer from this weakness, but in all the works cited above, extracting e.g. the probability 
density function (PDF) of the observed angular distance at a fixed redshift requires numerical 
simulations which, because of their computational cost, lack of flexibility. A practical solution 
was proposed with the sGL method of Kainulainen and Marra [83-85] , in which weak-lensing 
simulations have been maximally optimised so that generating 10^ mock observations only 
takes a few seconds. This method has been applied to forecast to which extent future 
SN observation campaigns, e.g. with the Large Synoptic Survey Telescope (LSST), would 
be able to constrain cosmological parameters from the moments of the distribution of SN 
magnitudes [86-89]. 

The goal of the present work is to propose an analytical and a priori non-perturbative 
framework for determining the statistical impact of small-scale structures on light propagation. 
Possible applications are the analysis of the bias and dispersion induced by these structures 
on cosmological observables, non only for distances measurements but also, e.g., cosmic shear. 
The main idea is that, on very small scales, the matter density held (i.e. the source of 
lensing) can be treated a white noise, giving to lensing a diffusive behaviour. The equations 
of geometric optics in curved spacetime then take the form of generalised Langevin equations, 
which come with the whole machinery of statistical physics. Indeed, similar approaches have 
been exploited in other domains of physics [90, 91], e.g., for describing the secular evolution 
of the Solar system. This systematic treatment of lensing as a stochastic process allows us to 

^The typical physical size of a supernova explosion is on the order of a hundred astronomical units, which 
fixes the typical maximum cross-sectional diameter of the associated light beam. On such scales, the distribution 
of matter in the Universe cannot be considered smooth. 
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derive Fokker-Planck-Kolmogorov (FPK) equations for the PDF of the lensing observables, 
such as the angular distance, on which we will particularly focus in this article. 

The benefits of this new approach are multiple. Its analytical character potentially 
provides a better physical understanding of small-scale lensing, together with avoiding to rely 
on heavy ray-tracing simulations. It must be considered complementary to cosmic lensing 
due to the large-scale structure, with which it is planned to be merged in the future, in order 
to design a consistent multiscale description of lensing. Similarly to Refs. [86-89], we have in 
mind applications to a better characterisation of the matter distribution within the Universe. 
These various applications lie beyond the scope of the present article, which however proposes, 
as starters: (i) an extension of the KDR approximation, and (ii) an analytical calculation of 
the variance of the angular distance in an Einstein-Straus Swiss-cheese model. 

The article is organised as follows. Section 2 provides a theoretical lensing toolkit, which 
contains all the necessary material exploited in the remainder of the article, in particular the 
Jacobi matrix and the optical scalars. Sections 3 and 4 are the heart of our approach: the 
former presents our fundamental hypotheses; the latter derives the FPK equations governing 
the PDF of the Jacobi matrix and of the optical scalars. Section 5 deduces general analytical 
results from the FPK equations, in particular regarding the hrst two moments of the PDF of 
the angular distance. In order to test our formalism, we apply it to a Swiss-cheese model, and 
confront the associated predictions to numerical ray-tracing results in Section 6. Section 7 is 
finally devoted to a second check of our calculations, based on the numerical integration of the 
Langevin equation using the stochastic Euler method. It sheds some light of the connection 
between the accuracy of our predictions and the Gaussianity of the sources of lensing. 

2 Propagation of narrow light beams: two complementary formalisms 

Consider a narrow light beam, that is an infinitesimal bundle of null geodesics, converging at 
an observation event O. Among the geodesics of the bundle, we arbitrarily pick a reference 
ray x^{v), where v is an affine parameter along the ray. The associated tangent vector 
= dx^/dv represents the wave four-vector of the light beam. If we choose k as past oriented 
(so V increases from O to the source), then the (cyclic) frequency measured by an observer 
crossing the beam with four-velocity u is uj = u^k^. In this article, we set by convention v = 0 
at O, and normalise all frequencies with respect to the observed one lOq = {u^k^)\o = 1- 

The behaviour of any ray x^{v) of the beam, relative to x^{v), is characterised by its 
connecting vector = x^ — x^. If an observer at x^{v) projects the beam on a screen, 
spanned by the Sachs basis (see Appendix A), then the relative position of the two light spots 
associated with and x^ is a Euclidean two-dimensional vector 

2.1 Jacobi matrix 

The first standard tool for describing the effects of gravitational lensing is the Jacobi matrix, 
whose evolution with light propagation is a seeond-order linear differential equation. 

2.1.1 Definition 

The Jacobi matrix is a 2 x 2 matrix T) = [Dab] which relates the physical separation 
(in screen space) between two rays with their angular separation .^^(0)—a dot denotes a 
derivative with respect to v —on the observer’s celestial sphere, according to 

= (2.1) 
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The determinant of X> thus represents the ratio between the beam’s cross-sectional area A{v) = 
at V with its observed angular aperture Oq = d^^. When evaluated at the source event 
{v = Vs), we recognise the definition of the (squared) angular diameter distance between the 
source and the observer 

detT>{vs) = ^ = Dl. (2.2) 

i Lq 

We recall that, if the number of photons is conserved during their travel from the source to the 
observer, then the angular diameter distance T>a is related to the luminosity distance—used 
e.g. in the Hubble diagram of SNe—by the distance duality relation 

Di^ = {1 + z)‘^Da, (2.3) 

which involves the redshift z = (wg — uJo)/oJo between the emitted and observed frequencies. 

The other three degrees of freedom of T) encode the deformations of the light beam, 
i.e. the deformations between the intrinsic source’s shape and the observed image. This 
information is conveniently extracted from X> by the decomposition given in Appendix A. 

2.1.2 Evolution: the Jacobi matrix equation 

Because describes the relative behaviour of two neighbouring light rays, its evolution with 
light propagation (i.e. with v) is inherited from the geodesic deviation equation; it results 
into the following second-order linear Jacobi matrix equation [92] 

ii) = 7l{v)'D{v) (2.4) 

where TZab = is called the optical tidal matrix, and (s(^)a=i ,2 denotes the 

Sachs basis. The optical tidal matrix is symmetric due to the symmetries of the Riemann 
tensor Rpupa- The decomposition of the latter into a Ricci (trace) part and a Weyl (trace-free) 
part implies, for the optical tidal matrix. 


7i = Agl2 + yV, (2.5) 

I 2 standing for the 2x2 unity matrix, while 

( 2 . 6 ) 

Was = C^,p„s^Xk'^kPs% (2.7) 

where R^^ and Cp^pa denote respectively the Ricci and Weyl tensors. It is straightforward to 
check that W is trace free, and can thus be written as 

w = , with Wi + m^w = -^Cp,p„is>; - i 4 )rF(s? - i^) (2.8) 

The Ricci term, on the one hand, is directly related to the local energy-momentum density via 
the Einstein equation, ^ = —ATrGTpuk^k^ < 0 (under the null energy condition); it translates 
the isotropic focusing effect caused by smooth matter enclosed by the light beam. The Weyl 
term, on the other hand, essentially encodes tidal distortion effects, due to matter outside the 
beam, which tends to shear and rotate it. 
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The initial conditions {v = 0) for Eq. (2.4) are by definition [see Eq. (2.1)] 


(2.9) 

( 2 . 10 ) 


X>(0) = O 2 

^(0) = l2, 

so that, near the observer (u —)• 0), the Jacobi matrix admits the expansion 


3 

T>{v) = vl2 + l.Tto + 0{v^). 


It also implies, using that for any matrix M, det(l + eiW) = 1 + etriW + 

,,3 


DAiv) = v + ^^o + Oiv^). 


( 2 . 11 ) 


( 2 . 12 ) 


2.2 Optical scalars 


A standard alternative to the Jacobi matrix consists in a set of optical scalars, describing the 
deformation rate of the beam rather than net transformations. The resulting light propagation 
equations (Sachs equations) are a set of first-order nonlinear equations. 

2.2.1 Definition 

The deformation rate of the light beam is naturally defined by a logarithmic derivative of the 
Jacobi matrix, namely through 

S = 'D'D \ (2.13) 

This deformation rate matrix can be shown to be symmetric, because of the symmetry of 72-, 
and is thus decomposed as 


S = 


9 0 
0 9 


+ 


— <Tl (T2 
(T2 (Tl 


(2.14) 


where 9 and cr = ui + \(T 2 are the optical scalars, respectively called the expansion rate and 
the shear rate. The first one is directly related to the increase rate of the angular diameter 
distance, since d(lndetX>)/du = tr«S, i.e. 


9 = 


Da 

Da 


(2.15) 


2.2.2 Evolution: the Sachs scalar equations 

Inserting the definition (2.13) into Eq. (2.4) yields the evolution equation for «S, 

S + S^ = TZ, 

from which the Sachs scalar equations follow: 

9 + 9 ^ + jfTj^ 

a + 29a = W. 

Using that 9 = Da!Da., the above equation yields the so-called foeusing theorem 

bA = {^-\a\^)DA, 


(2.16) 


(2.17) 

(2.18) 


(2.19) 
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where we see that, while Ricci lensing has a direct focusing effect which tends to rednce 
Weyl lensing has a similar but indirect effect, via the shear rate. 

The initial conditions for the optical scalars are nontrivial, because X> vanishes for 
n = 0, which implies that S must have a pole at the observation event. Precisely, the initial 
behaviour (2.11) of the Jacobi matrix yields 

<S(n) = [I 2 + 0(^2)] [v I 2 + 0{v^)\ = v-^U + 0{v), (2.20) 

and we conclude that the initial conditions {v —?■ 0) for the optical scalars are 

e{v) = - + 0{v), ( 2 . 21 ) 

V 

a{v) = 0{v). (2.22) 

Hence only the expansion rate has a pole at n = 0, while the shear rate is regular. 

3 Small-scale lensing as a diffusion process 

We now focus on the specific issue of lensing caused by the small-scale inhomogeneity of the 
Universe, i.e, down to scales where the matter distribution experienced by the light beam 
cannot be considered a continuons medinm, but rather by a multitude of mass clumps that all 
slightly distort it. This situation is analogons to the Brownian motion of a particle snspended 
in water, where a macroscopic—continuous-medium—description of the liquid is no longer 
sufficient, and must be replaced by a semi-microscopic approach in order to account for the 
collisions between the particle and water molecules. 

The approach developed in the present article is based on this analogy. Just like in the 
standard treatment of the Brownian motion, where particle-molecule collisions are modelled 
by a stochastic force, we propose to introduce stochastic terms in the lensing scalars W. 
The equations governing light propagation will thus take the form of Langevin equations. 

3.1 Fundamental hypotheses 

We split the Ricci and Weyl lensing scalars experienced by the light beam into a deterministic 
part representing their average, slowly varying behaviour, and a stochastic part modelling 
their rapid fluctuations; 


= (3.1) 

ir = (>r) + j/r, (3.2) 

where (...) is an ensemble average, and = {b'W) = 0. All these quantities are in principle 
functions of the affine parameter. Note that, despite the notation, bS^ and bW are not 
necessarily small with respect to {t%) and (i^), they are not dealt with as perturbations. The 
deterministic components can be thought of as the optical properties of an average universe, 
in the sense e.g. of Ref. [54]— a notion which may not coincide with a spatial average, or with 
a FL model. 

We now make the following hypotheses: 

Azimuthal symmetry about the beam. We suppose that the Universe is statistically 
homogeneous and isotropic, which implies statistical symmetry with respect to rotations 
about any light beam. This motivates us to assume that the direction along which a 
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beam is sheared is independent from the shear amplitude. It is also independent from 
Ricci focusing. In other words, decomposing the Weyl lensing scalar as W = 
we assume that j3 is statistically independent from \if\ and However, we emphasize 
that \'if\ is not independent from 

Statistical isotropy. We snppose that the Universe has no preferred (spatial) direction, 
which implies that /3 must be uniformly distribnted in [0, vr]. As a consequence, 

{W) = {\W\) = 0, (3.3) 

where we have also used our first hypothesis. We can thns omit the 6 in the stochastic 
part of W. Furthermore, for any v, w 

{6^{v)W{w)) = {St%{v)\W{w)\) = 0, (3.4) 

{Wi{v)W 2 {v)) = ^ (|#'(u)p) (sin4;0(u)) = 0. (3.5) 

White noises. Because they model rapidly fluctuating functions, the coherence scale of 
and W is mnch smaller than the typical evolntion scale of the Jacobi matrix, of 
the optical scalars, and than the typical distance between the sonrce and the observer. 
Therefore, they can be considered white noises, i.e. 5-correlated Gaussian random 
processes^ , with 


{6^{v)6^{w)) = C^{v)6{v — w) (3.7) 

{Wa{v)Wb{w)) = C^{v)6ABdiv - w), (3.8) 


where the 6ab in Eq- (3.8) comes from statistical isotropy. The functions C^, C-p- shall 
be called the covariance amplitudes of Ricci and Weyl lensing. Gaussianity, which is 
motivated by the central limit theorem, ensures that S^{v) [resp. W{v)] and 5t%{w ^ v) 
[resp. W{w ^ u)] are not only uncorrelated, bnt also independent. 


Physically speaking, the covariance amplitude Cx of the white noise X{t) modelling 
a physical process Xphys(t) must be understood as Cx (JXphys)^Atcoh, where JXphys is 
the typical fluctnation amplitnde of Xphys, while AGoh is the scale on which it remains 
coherent. For classical Brownian motion, this scale corresponds to the duration of a typical 
particle-molecule collision; in gravitational lensing, it will represent the typical extension of 
a gas cloud/dark matter halo (Ricci lensing), or the affine-parameter length over which the 
beam undergoes the tidal influence of a given deflector (Weyl lensing). 

In principle, the deterministic components {M) and {W) could also allow for the large- 
scale structure of the Universe (cosmic voids, walls, and filaments). For simplicity, we do 
not consider this possibility in the present paper, and focus our attention on the rapidly 


® A random process t \-A X(i) is Gaussian if any of its finite-dimensional probability distributions is a 
multivariate Gaussian, 


Pti,...t„(*!,.. .Sn) oc exp -- ^ XiCC 




(3.6) 


where Cij = C{ti,tj) = {X{ti)X{tj)) is the covariance of the process, and C~^ denotes its inverse. A white 
noise corresponds to the limit where C{ti,tj) oc Hence, for a white noise, X{ti) and X{t 2 yf ti) are 

independent. 
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fluctuating terms. It will be convenient, in the following, to gather them into a 3-dimensional 
noise vector N such that 

= (,5^,#i,#2)- (3.9) 

We also introduce the diffusion matrix Q of N, deflned by^ (^N(v)(w)) = Q{v)6{v — w), 
which here reads 

Q = diag(C'^, C^, C^). (3.10) 


3.2 Langevin equation for the Jacobi matrix 

The Jacobi matrix equation (2.4) reads 

T> = {^)T> + {d^ + W)T>, (3.11) 


where we have separated the deterministic and stochastic terms on the right-hand side. It is 
analogous to a system of coupled harmonic oscillators with fluctuating stiffness. Some further 
insights on this dynamical system can be obtained thanks to a Hamiltonian formulation 


T^ab = Vab = 
'Pab = — 


dH 
QVab 
5Hjac dH 


(3.12) 


Pab 


QVab 


+ J^ABiv) 


with 

H =^tT{V^V - , Ar={5m2 + W)V, (3.13) 

and where the Hamiltonian H encodes only the non-stochastic part of the process. Such a 
dynamics is very similar to the integrable systems with stochastic perturbations discussed 
e.g. in Ref. [90], except that (i) due to the explicit u-dependence of H, through {^), the 
unperturbed system is not fully integrable; and (ii) the stochastic term JV contains the 
variable X>: the noise is multiplicative. This analogy with dynamical systems in statistical 
mechanics also provides a nice interpretation of the deformation rate matrix «S: as a Ricatti 
variable associated with X>, it defines the so-called Kolmogorov-Sinai entropy of the random 
process, /iks = tr(«S). 

Let us now put the Jacobi matrix equation in the form of a first-order Langevin equation, 
which will be useful for deriving the associated Fokker-Planck-Kolmogorov equations in Sec. 4. 
For that purpose, we first need to vectorise the Jacobi matrix as 


D = (T>a)ae{i... 4 } with Vab = D 2 (a-i)+b\ (3.14) 

in other words, we represent the couples of matrix indices (AB) by one single index a, so 
that 1 = (11), 2 = (12), 3 = (21), 4 = (22). We then construct an 8-dimensional vector 
= {D, D), whose dynamics is described by the Langevin equation 


dJ 

du 


M{v)J{v) + LjUJ)N{v). 


(3.15) 


where the drift matrix is 

(3. 


'^Equivalently, the diffusion matrix can be defined from the increments of the Brownian motion B associated 
with N, i.e. such that AB = Ndv. Between vi and V 2 , the increment of B is AB = B{v 2 ) — B{vi), and its 
variance reads = QAv, with Av = V 2 — vi. 


M = 


O4 I4 

{^) I4 O4 


- 8 - 









and the noise-mixing matrix reads 


- 1 

o 

X 

CO 

_ 1 


1 - 

CO 

X 

0 

1 _ 

Di —Di Ds 


7^11 —Dll D 21 

D 2 —D 2 D 4 

= 

V 12 —Di2 D 22 

Dz Dz Di 


V 21 D 21 Dll 

D^ ZI 4 D2_ 


P 22 D 22 Di2_ 


Equation (3.15) is linear, with a mnltiplicative noise. 


(3.17) 


3.3 Langevin equation for the optical scalars 

A similar procedure can be achieved for the optical scalars. The Sachs equations (2.17-2.18), 
together with the relation (2.15) between the angular distance and the beam’s expansion rate, 
form the system 


Da = ODa, (3.18) 

9 =-0^+ (3.19) 

& = -2ea + W, (3.20) 


which, defining the 4-dimensional vector = {DA, 6 ,o'i,a 2 ), becomes the Sachs-Langevin 
equation 

1 c* 

— = F{v,S) + L,,,iN{v), (3.21) 


where the drift term reads = {9 Da, 
mixing matrix is 


L 


seal — 


- kl" + {^) 

0 0 o' 

1 0 0 
0 10 ’ 

0 0 1 


—29ai, —20(72), while the noise 


(3.22) 


Contrary to Eq. (3.15), Eq. (3.21) has a nonlinear drift term (which reflects the nonlinearity 
of the Sachs scalar equations), but its noise is additive, in the sense that the stochastic term 
Lscai^AT is independent of the variable S. 


4 The lensing Fokker-Planck-Kolmogorov equations 

The presence of stochastic terms in the optical equations gives a diffusive behaviour to the 
lensing observables, which can be quantified by their PDFs. When a dynamical system is ruled 
by a Langevin equation, its PDF in phase space satisfies a partial differential equation called 
the Fokker-Planck-Kolmogorov (FPK) equation. In § 4.1, we recall the general procedure to 
derive the FPK equation associated with a Langevin equation; we then apply it to the Jacobi 
matrix (§ 4.2) and to the optical scalars (§ 4.3). 

4.1 Prom Langevin to Fokker-Planck-Kolmogorov 

Consider the following general Langevin equation governing the evolution of a n-dimensional 
random process te^X{t), 


= f{X,t) + L{X,t)N{t), 
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dW 

dt 


(4.1) 








where the ra-dimensional vector f and the n x n matrix L are deterministic, while AT is a 
white noise. One can easily see that both our Langevin equations (3.15) and (3.21) have this 
form, the affine parameter playing the role of time t, and the random process being either J 
or S. The mathematical difficulty of Eq. (4.1) is that it cannot be treated with the ordinary 
theory of differential equations, because N{t) is discontinuous everywhere. In general, the 
solution of Eq. (4.1) is not unique, even for a given realization of N. 

A standard approach [93-98] consists in introducing the Ito calculus, the main properties 
of which we summarise below. One can formally integrate Eq. (4.1) as 


X{t)-X{to) 


f f{X,t) dt + 

Jto 


ho 


L{X,t) N{t)dt, 


(4.2) 


where the second integral requires particular attention, because the Riemann or Lebesgue 
definitions cannot apply, due to the unboundedness and discontinuity of N. First, it must be 
reformulated as a Stieltjes integral 


L{X,t) dB 


'to 


(4.3) 


where S is a Brownian motion, i.e. a stochastic process whose any increment = B{tk+i)~ 
B{tk) is a zero mean Gaussian random variable with variance = Q{tk,tk+i)Atk- 

Q is called the diffusion matrix of B. The white noise N is thus considered a formal derivative 
of the Brownian motion B, i.e. dB = Ndt. One possible definition for the integral (4.3) 
follows the so-called ltd stochastic prescription [99], 



L{X,t) dB 


n—1 

lim '^L[X{tk),tk] [B(4+i) - B(4)] 

n^oo 

k=0 


(4.4) 


This definition leads to some modifications with respect to ordinary differential calculus 
when B is involved. For example, it can be shown by calculating explicitly the Ito integral of 
BidBj that d{BiBj) = BidBj + BjdBi + Qijdt, which implies 

dBdB^ = Qdt. (4.5) 


The above quantity is thus of order 1 in dt, contrary to what we would naively expect by 
replacing dB by Ndt. Equation (4.5) is the most important rule of the Ito calculus. As a 
consequence, the first-order Taylor expansion of any function f{t,X) must actually include 
second-order terms oc dXi dXj , since 

dX = f{X,t)dt +L{X,t)dB, (4.6) 

contains dB. More precisely, 

(4.7) 

(4.8) 

which is known as the ltd formula [93-99] . 


11 I 


d‘^(b 

2 


' dt 2 dXidXj 


d‘^6 ^ \ , dd) 

LikQkiLji \ dt -|- dXi, 
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From the Ito formula, one can deduce the Fokker-Planck-Kolmogorov (FPK) equation 
governing the PDF p{t‘,X) of the stochastic process X(t). The derivation [94, 96] relies on a 
trick which consists in inserting Eq. (4.8) in the time derivative of the expectation valne of an 
arbitrary function 


{cP)it) = j (t>{t,X)p{t-X)d^X, (4.9) 

which, after a few integration by parts, yields 

= - ay + 5 aiax, {lHt:x)Qm^(x.t)],^ p(t; x)}. (4.io) 

The first term on the right-hand side is a drift term, it drives the global displacement of the 
probability packet, while the second is a diffusion term, which tends to spread it. With this 
summary of textbook results [93-99] we wish to emphasize that the derivation of the FPK 
equation requires the noise to be white, i.e. N = dB/dt where S is a Brownian motion, so 
that the ltd calculus can be applied. The hypotheses formnlated in § 3.1 are therefore crucial 
for this formalism to be applicable. 


4.2 FPK equation for the Jacobi matrix 

Let us now derive the FPK eqnation governing the PDF of the Jacobi matrix. Applying the 
general formula (4.10) to the Langevin equation (3.15) leads to the following equation for the 
PDF p{v, J), 

dv ^ ~Wa + 2 dJadJk 

where the indices a, 6 rnn from 1 to 8. Using the explicit expression (3.17) of Tjac, we can 
write the matrix involved in the diffusion term as 


Lja,cQL 


T 

Jac 


O 4 O 4 

04 r ’ 


(4.12) 


where the components of the 4x4 symmetric matrix F are 

Til = {C^ J- + C^V^i 

ri2 = {C^ J- C^)ViiVi2 + Cy^V2lD22 = F 21 

Tis = C^V2iDii = F 31 

Fi4 = {C^ — C^)ViiV 22 + C'#'D2iT’i2 = F41 

r22 = (C^ + Cw)'^‘l2 “k C' 5^'^22 

r23 = {C^ — Cy^)Vi2T>2l + C^ViiD22 = F 32 
F 24 = C^Vi2'B22 = r42 
Tss = + C^)'D2i + 

F34 = {C^ J- 0^)1)211)22 + Cp'T)iiT)i2 = F43 
F44 = {C^ J- C^)T)22 + C^'D\2- 


(4.13) 
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A few calculations and reorganizations yield the following explicit form of the FPK equation 
ofp(u; J) =p{v,T>,T>), 


dp 

dv 


—i'AB 


dp 

dVAB 


{M) Vab 


dp 

dT>AB 


+ 2 ^AE^CF + C^{6aC^EF — £AC£Ef)] TdEBT^FD 


d'^p 

dVAEdVcD 


(4.14) 


where eab is the two-dimensional antisymmetric matrix with £12 = 1. Equation (4.14) can 
also be rewritten in an elegant formal way as 


dp 

dv 


— tr 



tr I X> 







C-p- det 



p (4.15) 


which involves in particular the 2x2 matrix differential operator 


—~ I = T^ca —=-• 

dVjAB dVcB 

Finally, the boundary condition for Eq. (4.14) is deduced from the initial conditions 
(2.10), and reads 



(4.16) 

(2.9), 


p(0;X>,X>) = ,5(X>)5(X)-l2). 


(4.17) 


4.3 FPK for the optical scalars 

Regarding optical scalars, starting from the Langevin equation (3.21), one can derive the 
following FPK equation for p{v; S) = p{v] Da, 0, cti, CJ 2 ), 


dp 

dv 


dFaP 1 d'^ 

dSa ^ 2 dSadSfi 


(-hscal SRgcal) ap P 


(4.18) 


where a, (5 run from 1 to 4, and where the diffusion term reads 


LscalQL 


T _ 
seal 


0 0 0 0 
0 0 0 
0 0 0 
0 0 0 


(4.19) 


It follows that Eq. (4.18) takes the explicit form 


dp ^dDxp d 

Ih “ aoT di 


’(«" + kl 


P 




d(T2p \ 
d(J2 ) 


+ 


Cs^d'^p CiA / d‘^p d'^p 
o ' o I ^^ 2 ' 


(4.20) 


The initial condition for 6 being singular, it is not possible to write a boundary condition for 
Eq. (4.20) as we did for Eq. (4.14). 
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5 General analytical results 


Because it is a partial differential equation, the FPK equation is generally impossible to 
solve analytically, except in a few known special cases [98]. Nevertheless, it can be used to 
derive evolution equations for the moments of the PDF, some of which are solvable. In this 
section, we derive some general analytical formulae on the moments of lensing observables. 
The results for the Jacobi matrix (§ 5.1) and for the optical scalars (§ 5.2) will turn out to be 
complementary, and used for deriving an evolution equation for the variance of the angular 
diameter distance in § 5.3. 

5.1 Moments of the Jacobi matrix distribution 

The Jacobi matrix formalism has this considerable advantage on the optical scalar formalism 
that it enjoys a linear Langevin equation. Despite the fact that its noise is multiplicative, 
this implies that all the moments of order-n of the PDF of satisfy a closed system of 
differential equations. It is not the case when nonlinearities are present, in which case emerges 
a hierarchy of equations, where the evolution of the lower-order moments depends on moments 
of higher-order. 

5.1.1 Order-one moments 

Let us start by deriving the evolution equations for the expectation values (X>) and {T>). We 
proceed by multiplying the FPK equation (4.14) by Vjj (or D/j) and then integrating it with 
respect to T> and T>. For Vjj, this procedure yields 


^ I Vjjpd^T>d^T> = - I d^Vd^T) - {M) I d^Vd^t) 

1 f 

+ o / —{[Cgg dAEdcF + {5 ac^ef — £ac£ef)]T^eb"Dfdp} d^T>d^T>. 

2 J dVABdVcD 

(5.1) 


The left-hand side is clearly d{'Dij) /dv. On the right-hand side, the hrst term can be 
integrated by parts to give (T>/j); the other two vanish since they can both be written as the 
integral of a derivative with respect to Vab- Equation (5.1) is thus simply 


d(^) 

du 


(P) 


as one can intuitively expect. 

The same method applied to Vij leads to 


dp) 

dv 


= m p) 


so that the expectation value of the Jacobi matrix reads 


d^ (X>) 
du^ 


m P) ■ 


(5.2) 


(5.3) 


(5.4) 


Note that this result could also have been obtained by directly averaging the Sachs-Langevin 
equation. However, this naive method would not work for higher-order moments, which is 
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why we preferred to directly use a rigorous technique for deriving the evolution equation for 
the expectation value of X>. 

It is tempting to conclude that the average angular diameter distance (Ha) satisfies 
Eq. (5.4) as well, but such an assertion would be wrong, because Da = VdefD is a nonlinear 
function of the components of the Jacobi matrix. 

5.1.2 Order-two moments 

We apply the same method to get evolution equations for the order-two moments of X>. This 
leads to the following closed system of equations 

{DabDcd) = {DabDcd) + {DabDcd) (5.5) 

du 

-^{DabDcd) = {DabDcd) + (■^) {DabDcd) (5.6) 

du 

■^{DabDcd) = (■^) (jPabDcd) + {Dab'Dcd)'^ + {DabDcd) 

-I- {SacSef — ^ac^ef) (DebDfd) (5.7) 

which consists of 10 + 16 + 10 = 36 independent equations for the quantities (DabDcd), 
{T>abDcd) and (VabDcd)■ By combining the second derivative of Eq. (5.5) with the 
derivative of Eq. (5.6) and Eq. (5.7), we can eliminate the moments {DabDcd) and {DabDcd), 
in order to end up with a closed system for {VabDcd), 

{DabDcd) = 4 {^) ^ (DabDcd) + 2 ^^ {DabDcd) 

du'^ du \ av J 

+ 2(7;^ {^AC^EF — sac^ef) {DebDfd) , (5.8) 

which consists of 10 independent third-order differential equations. We shall not try to solve 
this system, but rather extract from it information on the angular distance. 

5.1.3 Application to the squared angular distance 

The square of the angular distance is the determinant of X>, hence quadratic in its components. 
Its expectation value. 


(Hi) = (detX>) = (H 11 H 22 ) - (H 12 H 21 ), (5.9) 


is therefore ruled by Eq. (5.8). Applying it for ABCD = 1122 and ABCD = 1221, we have 


(H11H22) = 4 (^) — (H11H22) + 2 (H11H22) — (-^a) ) 

^ (H12H21) = 4 ^ (H12H21) + 2 (H12H21) + 2 C^ (Hi) , 

which, by subtraction, yields the following equation for {D'j^ only, 



(5.10) 

(5.11) 


(5.12) 
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To our knowledge, it is the first time that such a general exact equation for the evolution of 
the dispersion of the angular distance in an inhomogeneous universe is derived. 

Solving this differential equation requires initial conditions for (T^a) 
second derivatives. They are easily obtained from the Taylor expansion (2.12) of Ha for 
u —)■ 0, 

(Hi>(0) = 0, 1^(0) = 0, ^!i^(0) = 2. (5.13) 

Equation (5.12) can also be elegantly rewritten in terms of a variable x defined by 


Ax = 



(5.14) 


where Dq{v) is the background angular distance, i.e. satisfying Dq = Dq. It is indeed 
straightforward to show that the differential operator involved in Eq. (5.12) reads 




, d dm 

du du 


= H 


-4 




a 


-2 


(5.15) 


so that 

^ (^) = 2^0 - 2C';r) {Dl ). (5.16) 

Though formally simpler, this alternative form of Eq. (5.12) cannot be used for numerical 
integration, because x is singular at the observation event— Do{vo) = 0—usually chosen as 
initial condition. 


5.1.4 Expectation value of a general function 

More generally, by multiplying the EPK equation with an arbitrary function F(T>,T>) and 
integrating the right-hand side by parts, we obtain 


d{F) 

du 


^ dF \ 


Fab 


dF \ 
dFAB / 


1 / F \ 

+ O ( [Csi ^AE^CF + {SaC^EF — £AC^Ef)]FebFfD —T TTr ) • (5.17) 

2 \ dVABOFcD / 


If F is an order-ri monomial of the form F = with p + q = n, and where stands 

for any product of p components of the Jacobi matrix, then the left-hand side of Eq. (5.17) 
is d(P^P'?)/du, while the three terms on the right-hand side are respectively of the form 
(pP-ipPi), (pP+ip9-i)^ and (pP+2p<?-2), so they are all order-n moments. This confirms 
what we claimed in the introduction of this section, namely that order-n moments form a 
closed system of differential equations. 


5.2 Moments of the optical-scalar distribution 

Contrary to the Jacobi matrix, the optical scalars satisfy a nonlinear Langevin equation. An 
important consequence on the associated EPK equation (4.20) is that it generates an infinite 
hierarchy of evolution equations for the moments of the distribution p(v; S). For instance, if 
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one is interested in computing the average angular distance {D\), then Eq. (4.20) generates 
(using the same technique as in § 5.1) 


- (ZIa) = {9Da) 

^{ODa) = -{ kl' Da) + m {Da) 
k|2 Da) = -3( k|2 BDa) + 2Cy^ {Da) , 


(5.18) 

(5.19) 

(5.20) 


where the evolution of an order-n moment systematically involves order-(n -|- 1) moments. 
Clearly, such a system cannot be solved analytically, and reqnires a perturbative approach to 
be dealt with. A first possibility consists postulating a closure relation for the hierarchy at a 
given order, but such a method does not seem particularly adapted to the present situation, 
becanse the physical meaning of the underlying approximation is unclear, and therefore poorly 
controlled. 

We choose instead to perform a perturbative expansion with respect to the shear rate cr, 
that we assume to be a small quantity. In the following, we focus on the average angular 
diameter distance {Da), and determine its evolution at first and second order in |crp. 


5.2.1 First-order perturbative expansion 

We decompose the angular distance and the expansion scalar as 

Da = Dq + Di, 

9 = do + 01, 


(5.21) 

(5.22) 


where Dq, already introduced in § 5.1.3, is the solution of Dq = {M) Dq, and 9o = Dq/Dq 
is the corresponding expansion rate; both are deterministic quantities. We assume that the 
stochastic quantities Di,0i are small, in the sense that their probability distributions are 
concentrated on values much smaller than Do,9q respectively. 

We then expand Eq. (5.18) and the following two equations, generated by FPK, 


^^{9) = -{9^)-{\a?) + 


(k|2> = -4(0k|k +2C1^, 


at first order in Di,6i,\a\ 


dn 

which gives 

d^ 

du 

d(0i) 


dv 

d(k|^ 

dv 


= 9o {Di) + {Oi) Do 
= -200 {9i)-{\a\^) 


= -40o(k|^) -I- 2C^, 


whence 


r d,.i 

Df, Jo r*o(^i) 


nvi 


dv2 


rv2 


dus Dq{v3)C^{v3) < 0, 


(5.23) 

(5.24) 

(5.25) 

(5.26) 

(5.27) 

(5.28) 
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which represents the relative correction between {D\) and Dq at first order in Weyl lensing. 
Note that, again, the above result naturally exhibits the integration measure D^^dv = dx, it 
can therefore be rewritten as 


dx^ 


-2DlC^. 


(5.29) 


5.2.2 The shear rate at first 

From Eq. (5.27) we have deduced 

(kP) 


order 

the following expression for the variance of the shear rate, 
4 


= 2 / drc 
Jo 


Dojw) 

Do{v) 


Cw{w), 


(5.30) 


at first order. Let us simply mention that, to this order of approximation, we can easily obtain 
the full PDF of a. Linearizing the second scalar Sachs equation (2.18), we indeed get 


d = -lOoa + W + 0{a^). 


(5.31) 


which is identical to the historical Langevin equation for diffusion. The associated FPK 
equation for the PDF Pa{v, a) is easily shown to be 


dpa _r.a ( dcriPa da2Pa \ CV ( d‘^Pa 
dv ° \ dai da 2 ) 2 da 2 ^ J 


(5.32) 


It can be solved by (i) using a polar description for a = ai + io '2 = |o'|e^'^, then (ii) using 
the statistical isotropy assumption that implies Po-(u; cii; CT 2 ) = f{v,\a\), and (hi) performing 
simple changes of variable to recover a standard diffusion equation. The result is a Gaussian 
distribution, describing a 2-dimensional random walk with nonconstant diffusion coefficient. 


Pa{v;a) 



1^1' ] 


(5.33) 


where (|cr|^) is given by Eq. (5.30). 

5.2.3 Second-order perturbative expansion 

In § 5.1.3 we derived an evolution equation for (D^), while § 5.2.1 provided an expression for 
(Da). Subtracting the results should therefore lead to the variance of the angular diameter 
distance. However, the first-order expansion performed in the previous paragraphs is not 
sufficient for that purpose. This can be understood the following way: if D\ = Dq -|- 5D, then 

var(DA) = (Dl) - {Da? = {SD? - {SD? (5.34) 

involves second-order quantities, neglected in § 5.2.1. In this paragraph, we therefore expand 
the equations governing the evolution of (Da) up to second order in (|(t|^), i.e. formally up 
to second order in C^. 

We start back from Eqs. (5.18), (5.19), (5.20), which can be gathered as 

S (^) = ^^a) + 3D5( lap Da{9 - 9q)). (5.35) 
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The difficulty now consists in evaluating the last term. First note that, since it is already a 
second-order quantity, 



(lap Z4 a( 0 - 0o)) = Do{ |ap 0i) + 0(C|.). 

(5.36) 

Let us then write 

(|ap0O = (kP)<0i> + (kp0>-(kP>(0>. 

'- V -" 

(5.37) 

— UCT 

The first term on the right-hand side can be expressed using the hrst-order results of 
which, using the x variable, take the simple form 

§ 5.2.1, 


{\-n = -Do^^ + OiCl.), 

(5.38) 


(«.> = Oo"^+0(C|,). 

(5.39) 


Evaluating the cross-correlation term can be achieved by using again the hierarchy 
of moments generated by the FPK equation. Combining Eqs. (5.23), (5.24) with 


|a|2 ) = -5(02 |a|2 ) - ( |a|^ ) + {^){ |ct|2 ) + 2C^{9), (5.40) 

we get 

f,, + moTe. = { lup >' - ( |ct|" ) + 0{C^), (5.41) 

where we have expanded the higher-order correlator (0^ |(t| ) as 0 q(|(t| ) -|-20o(0i \a\ ) +0{C^). 
Now, by comparing the evolution equations for and (|cr|^), which are 

1*^1^ 1^1^ )( I'^l^ ) + Wf )> (5-42) 

Wf ) = -8(^ Wf ) + 8C^{ \a\^ ), (5.43) 

we conclude that (Icrl^) = 2{\a\^)‘^ at leading order. Note that this result coincides with 
the predictions of the Gaussian distribution (5.33) obtained for a in the previous paragraph. 
Hence Eq. (5.41) is solved as 


F,, = -Dq-s r du; Dl{ lap >' + ©(C^) 
Jo 


= -D 


-6 

0 



+ 0{C^^), 


(5.44) 

(5.45) 


where we used Eq. (5.38). The lower bound “o” of the latter integral is formal, because 
variable x is singular for u = 0. This was the last missing piece to the differential equation 
governing the evolution of (T>a) at second order in 


({Da)\ 

dx3 \ Do J 


+ 2Cy^Dl 


(Da) 

Do 


-3- 


^ '^Da 
dx dx^ 



dx^ 


2 


(5.46) 
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In terms of an expansion of the form Da = Dq + Di + D 2 , and defining the second-order 
mean correction = (D 2 ) /Dq to the angular distance, the above result reads 




o Dp^ Da 

dx dx^ 



dx^ 


2 

< 0 . 


(5.47) 


5.3 Variance of the angular distance 

We now have enough material to propose an approximate evolution equation for the variance 
of the angular distance. On the one hand, we have obtained in § 5.1 the following exact 
equation for (D^), 

(5.48) 

On the other hand, the second-order Eq. (5.46) is easily turned into an equation for (Da) , 


d^ ({Da) 
dx^ \ Do 





dx^ 


2 

+ C>{C^). 


(5.49) 


By subtraction, we finally obtain 



(5.50) 

where we recall that dx = D^'^dv, and that the third derivative d^/dx^ is given by Eq. (5.15). 
We see that both Ricci lensing and Weyl lensing drive the variance of Da. This can be 
easily understood from the focusing theorem (2.19), where M is the main driving term, which 
explains why appears directly on the right-hand side of (5.50); W, on the other hand, 
affects Da only indirectly, via \(t\^■ It is the reason why d^d^^/dx^ oc (|(t|^) is also present on 
the right-hand side of Eq. (5.50). 

It is remarkable that this result on the variance of Da required the use of both the 
Jacobi matrix and the optical scalars. Although they are completely equivalent formulations, 
it would have been much more painful to derive Eq. (5.50) by using exclusively one of them. 

6 Application to a Swiss-cheese model 

The stochastic lensing formalism developed throughout Secs. 3, 4, and 5 depends on three 
free functions: the average Ricci focusing {M) (x), and the two covariances amplitudes C^{v), 
C-^{v) which need to be specified, or deduced from a spacetime model, in order to draw 
any physical conclusion. In this section, we propose an application of this formalism to 
Swiss-cheese (SC) cosmological models. Our goal is twofold: on the one hand, it provides an 
explicit example about how stochastic lensing can be applied, and of the involved calculations; 
on the other hand, it allows us to test its validity, by comparing its analytical predictions 
with the numerical results of a ray-tracing code for SC models, which was developed by one 
of the authors and used in Refs. [71, 100]. As a byproduct, we also obtain an improvement of 
the Kantowski-Dyer-Roeder approximation, which allows for shear. 
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6.1 The Einstein-Straus Swiss-cheese model 

We consider here an Einstein-Straus [23, 24, 101-103] SC model, where individual masses, 
whose vicinity is characterised by the Schwarzschild solution (or the Kottler solution, for 
a nonvanishing cosmological constant), are embedded in an expanding homogeneous and 
isotropic Universe, forming spherical holes within the Friedmannian cheese. This model aims 
at describing static, gravitationally bound objects, such as stars, galaxies, or clusters of 
galaxies, and is therefore more adapted to the problematic of small-scale inhomogeneities 
tackled here than LTB [104, 105] or Szekeres [105, 106] Swiss-cheese models. 


6.1.1 Spacetime geometry 

Let us briefly summarise the main geometrical properties of the Einstein-Straus model—more 
detailed explanations can be found, e.g., in our previous works [71, 100]. Consider one hole 
of the SC, whose centre is taken to be the origin of the coordinate system, without loss of 
generality. On the one hand, the metric of the exterior region is 


ds2 


-dT^ + a^{T) 


■ di?2 
1 - KR‘^ 


+ dU^ 


( 6 . 1 ) 


with dO^ = + s\v? 9dif^, K = cst, and where the evolution of the scale factor a with 

cosmic time T is ruled by the Friedmann equations, in particular 


/1 da \ ^ SttGpo / ao \ 3 K A 
\ a dT) 3 Va/ a2~*~3’ 


( 6 . 2 ) 


where po is today’s mean density of matter, modelled by a pressureless fluid. The cosmological 
parameters quantifying the relative importance of matter, spatial curvature, and cosmological 
constant in the expansion dynamics are respectively = 87rGpo/{3H^), = -KjiaHf, 

and Oa = A/(3i7^). The interior geometry is, on the other hand, given by the Kottler (or 
Schwarzschild-de Sitter) metric 

A 2 

ds^ = — A{r) dt^ + A~^(r) dr^ + dVt^ with A(r) = 1 — —-(6.3) 

r 3 

and where rg = 2 GM is the Schwarzschild radius associated with the mass M at the centre 
of the hole. 

The metrics (6.1) and (6.3) are glued together on a spacelike hypersurface corresponding 
a comoving sphere (the boundary of the hole), hence dehned by i? = = cst in terms 

of exterior coordinates, and r = rh(t) in terms of interior coordinates. The Darmois-Israel 
junction conditions [107-109] then impose 


rh(t) = o(r)i2h, (6.4) 

47r Q 

M = (6.5) 

Equation (6.5) must be understood as follows: the mass M at the centre of the hole is identical 
to the one that should be contained in the sphere of comoving radius R^, if the latter were 
homogeneously filled with the same comoving density po as the exterior. 
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6.1.2 Optical properties of each region 


Within the cheese, since the FL metric is conformally flat, light rays follow straight lines 
in terms of a suitable coordinate system. The cyclic frequency of the associated wave, as 
measured by a comoving observer (with four-velocity Ot), and normalised by the observed 
frequency at O, reads 


1 I I.T 


( 6 . 6 ) 


from which follows the relation between redshift z and affine parameter v for a light ray 
propagating through the cheese only. 


du _ 1 

dz ff (z) (1 -h z)^ ' 

Besides, the Ricci and Weyl lensing scalars are shown to be 

= -47rGuj^p(T) 
#FL = 0. 


(6.7) 


( 6 . 8 ) 

(6.9) 


Inside the hole, a light ray propagating in the 0 = 7r/2-plane admits two constants of 
motion, E = A(r)k^ and L = respectively associated with the stationarity and spherical 
symmetry of the metric. Their ratio defines the impact parameter b = L/E^ roughly equal to 
the closest approach radius r min ~ 6 of the photon trajectory if b ^ rg. The Ricci and Weyl 
lensing scalars read, in this case, 


— 0 






( 6 . 10 ) 

( 6 . 11 ) 


where /3 is the impact angle, corresponding to the angle between the plane of the trajectory 
and the hrst vector of the Sachs basis, as represented on Fig. 1. 


6.2 Effective optical properties 

Because of the intrinsically discrete nature of the SC model, we need to design an effective 
approach to be able to use the formalism developed in this paper. 

6.2.1 The Kantowski-Dyer-Roeder approximation 

The hrst set of effective optical properties for SC models was proposed by Kantowski [65] in 
1969, assuming that the mass clumps modelled by the central mass of the holes are extended 
and opaque, i.e., imposing a cutoff for the impact parameter b > 6minj which corresponds to the 
physical radius Tphys of the clump. This work was generalised in 1974 by Dyer and Roeder [69] 
in order to include the cosmological constant. The resulting behaviour at lowest order, that we 
shall call the Kantowski-Dyer-Roeder (KDR) approximation, can be summarised as follows: 

KDRl The relation between affine parameter v and redshift z is not signihcantly affected by 
the holes, so that Eq. (6.7) can still be applied in a SC model. 

KDR2 The effect of the shear, due to Weyl lensing in the holes, on the angular distance is 
negligible. In other words, i^KDR = 0. 
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Figure 1: Impact parameters in a Kottler hole. The grey disk is the intersection between 
the hole and the plane orthogonal to the wave-vector k at minimal approach, also spanned 
by the Sachs basis (si,S 2 ) there. The impact parameter b = L/E is approximately the 
minimal approach radial coordinate of the photon, and /? is the angle between the plane of 
the trajectory and the plane spanned by fc, si at minimal approach. It is also the angle 
corresponding the basis change which diagonalises the optical tidal matrix 7 ?.k in the hole. 


KDR3 Ricci lensing is the same as in the cheese but reduced by a factor a G [0, 1], called 
smoothness parameter, so that .^kdr = = —47rGa;^d/3(T). 


A detailed analysis of this approximation was presented recently in Ref. [71]. Hypoth¬ 
esis KDRl turns out to be valid up to terms on the order of the ratio rs/^h between the 
Schwarzschild radius of the central mass and the radius of the hole, which is very small in 
practice. Therefore, we will adopt KDRl for the remainder of this article. The relevance of 
KDR3 can be understood as follows: consider an interval [vmVn+i] of the light path, where 
Vn corresponds to the entrance into the hole number n, and Vn+i = Vn + to the entrance 
into the next one; the effective Ricci focusing over this interval can be defined as 


\ fVn+l 

= -- / M dv 

Au„ 




■^FL, 


( 6 . 12 ) 


where AuJj'^ is the fraction of light path spent into the FL region (between the exit from the 
hole n and the entrance into the hole n -|- 1) over which .^fl can be considered constant. This 
dehnes a local smoothness parameter = Au]j^/Au„. Interpolating the sequence (a„) on 
the whole light path yields a function a{v) which, after averaging over many lines of sights, 
dehnes a{v). 

In terms of the stochastic lensing formalism, we can thus identify 


(.^) = ^kdr- (6.13) 

As a consequence, the angular diameter distance predicted by the KDR approximation 
corresponds to Dq introduced in § 5.1.3, i.e. satisfying Dq = Dq. 


6.2.2 Effective Weyl lensing in a hole 

Numerical ray-tracing simulations in SC models [71] show that, while the KDR approximation 
satisfactorily reproduces the true Dx{v) relation for most lines of sights, some exhibit signihcant 
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deviations. Such discrepancies are due to Weyl lensing, neglected in the KDR approach 
(KDR2), but which we would like to include in the stochastic approach. It will be convenient, 
for that purpose, to first derive an effective expression for the Weyl lensing scalar in a 
single hole, defined as 

1 rvout 

Wes= - / #k(u) du, (6.14) 

^out '^in Jvin 

where Uin, Uout respectively denote the affine parameter at entrance and exit. 

Like in the KDR approach, we assume from now on that the central mass is an extended 
opaque object, whose physical radius rphys rs is thus a lower cutoff of impact parameters. 
As shown in Ref. [100], the radial coordinate r{v) of a photon propagating through the hole 
with an impact parameter b reads, at lowest order in rs/b, 

r{v) ^ ^/b‘^ + E^{v-v^y, (6.15) 


where Um denotes the affine parameter at minimal approach, and E k, Win ~ Wout- Moreover, 
if we neglect the growth of the hole between the photon entrance and exit, then 


l^out ~ Um Uin ~ 



Calculating the integral of Eq. (6.14) thus yields 




gme"^ 



2 

62rh 


e 


-2i/3 






e 


-2i/3 


(6.16) 


(6.17) 

(6.18) 


We see that, for b <C rh, the ratio between Weyl and Ricci lensing can actually be very large, 
oc (^h/&)^- It is the randomization of /3 which, in practice, drastically reduces the 
net impact of Weyl lensing on the angular distance. 


6.3 Calculation of the covariance amplitudes 

We now turn to the calculation of the statistical quantities C^, of the white noises which 
best reproduce lensing in a Swiss-cheese model. 


6.3.1 Statistical setup 

The randomness of our SC model is constructed in a way that—as originally formulated 
by Ref. [72]— “each ray creates its own Universe”. One realization of the various stochastic 
processes at stake thus corresponds to the disposition of successive holes on a photon’s 
trajectory, with random sizes, impact parameters, and separations. Expectation values (...) 
will be considered with respect to such realizations. As in Ref. [71], we make the following 
assumptions: 

• The properties (mass, size, impact parameters) of two different holes are independent, 
as well as the separation between different successive holes. 

• All the impact positions, within a given hole cross-section, are equiprobable. In other 
words, the impact angle {3 is uniformly distributed in [0,27r], and the PDE of the 
comoving areal impact parameter B is 

p{B) dB = \Rc<B < Rh] pf *^^2 ’ 
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where the squared bracket is 1 if the assertion inside is true, 0 if not; Rc denotes the 
comoving areal radius of the central matter clump, and i?h the comoving areal radius of 
the hole. We assume that the matter clump is static, i.e., its physical radius = aRc is 
constant, hence iic oc is not, contrary to i?h- 

• The distributions of both R\^ and are governed by the specific matter clumps that one 
wishes to model. For most of our theoretical results, they do not need to be explicitly 
specified. For numerical illustrations, we consider galaxylike clumps which all have the 
same physical density pc = 3M/(Airr^) = 3.47 x 10“^^kg/m^—this fixes the relation 
between Tc and M (hence i2h)—j and whose mass function is inspired from Ref. [110], 

p(M)dM oc (- ^ 5 ,, dM. (6.20) 


• The PDF of the comoving separation AypL between two successive holes is taken to be 
uniform, between 0 and 2 (Axfl), with 

(A»L> = 5^(Rk>. (6.21) 

3 1 — a 

This choice ensures that the mean smoothness parameter {a) = (Aufl/Au) is indeed a. 

6.3.2 Ricci-lensing covariance 

In reality, the Ricci and Weyl lensing scalars in a random Swiss-cheese model are not white 
noises: they have a self-correlation length on the order of the hole sizes. We here aim at 
determining the properties of the white noises which best reproduce the actual behaviour of 
and W. In the case of the Ricci covariance amplitude, this can be achieved by integrating 
Eq. (3.7) with respect to w, 


cMv) 


dw (6^(v)d^(w)} 


J dw {6^es{v)S^es{w)) , 


( 6 . 22 ) 


with 

= -^eff — ~ (6.23) 

and 5a{v) = a{v) — d. As mentioned above, the expectation value (...) is identified with an 
average over all possible realizations (r) of the SC, that is over the position, size, and impact 
parameter of each hole that is crossed by the light beam, 

1 ^ 

= Ymi ^ ^<5^efr('i')(^^eff(ti^)- (6.24) 

N^OO iV ^ ^ 
r=l 

For each realization (r), the complete light path through the SC can be split into 
elementary intervals R = [u„,u„+i], of affine parameter length Avn where, as before, Vn 
corresponds to the entrance into the nth hole. Within each interval, = 5^n is considered 
constant, and 5^n is independent of 6^rn if n ^ m. Hence, if we call 7(r)(^^) the elementary 
interval of (r) such that v G /(i.)(u), then there are two categories of realizations: those where 
w G as well; and those where w 0 /(j.)(u). The net contribution of the second category 

to the sum of Eq. (6.24) vanishes. 
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In order to calculate this sum, it is convenient to sort the realizations (r) in terms of 
the properties of The affine-parameter length Av of any elementary interval I can be 

decomposed into its FL and hole contributions as 

Av = Aufl + Auh = ^ ^Axh + , (6.25) 

where we neglected the global beam deflection in the hole part, and used the FL relation 
between affine parameter and comoving distance, even in the hole^. Av thus depends on the 
random parameters Axfl, Rin and B, which we regroup in a triple 11 = (AxfLj ^h, i?)- We 
now organise the sum of Eq. (6.24) in terms of the parameters 11 characterizing the interval 
containing v, which yields 

{6^esiv)^^es{'w)) = [ dll p(n|u G In) Prob(t(; G Iu\v G /n,n) (6.26) 


In the above equation, p(Jl\v G /n)dn represents the (conditional) probability that the 
interval In containing v has its parameters within dll around 11. It can be rewritten thanks 
to the Bayes formula as 


p{U\v G In) 


Prob(u G ln|n) 
Prob(u G I) 


X p(n). 


(6.27) 


where p(n) is the unconstrained PDF of 11, i.e. as provided by the assumptions of § 6.3.1. 
Simple geometric arguments show that the probability that v belongs to a given interval Ini 
with affine-parameter length At>(n), is 


Prob(u G 7n|n) oc Av, 


(6.28) 


so that the normalization factor in the denominator of Eq. (6.27) is simply Prob(u G /) oc 
(Au)jj, where the average is performed with respect to p(n). 

The second term in the integral of Eq. (6.26) represents the probability that w belongs 
to the interval In, given its parameters 11 and the fact that v already belongs to it. Again, 
simple geometry yields 


Prob(t(; G In\v G In, n) = ( 1 — 


|u — w\ 
Av 


Q{Av — \v — rcl), 


(6.29) 


where 0 denotes the Heaviside function. Gathering all the results, and using the expression 
of SMeS, we obtain 


{5^es{v)5MeE{w)) = {A-KCpobj^f [ dn p(n) —— 0(Au - ju - m;|) 

Performing the integration, plus the one with respect to w, finally yields 



/ Aufl 
V Au “ 

(6.30) 


(6.31) 


^This operation is justified by KDRl, which is very accurately satisfied in a SC model 
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in terms of the usual cosmological quantities. In the above equations, angle brackets denote 
averaging with respect to the mass function of the matter clumps, which rules the size of the 
hole they belong to via Eq. (6.5). Note that we get = 0 in both limits d = 0,1. This 
was indeed expected: for d = 0 the Swiss cheese is completely filled by holes, so that M = Q 
everywhere; for d = 1, we recover the strictly homogeneous FL spacetime, in which M 
everywhere. In both cases the fluctuation 8M vanishes. 


6.3.3 Weyl-lensing covariance 

Just like in the Ricci case, the covariance amplitude of the white noise which best 
reproduces Weyl lensing in a SC model is 

Cw{v) = \ ldw {W{v)W*{w)) ^'^Jdw (^\Wes{vWes{w)\ , (6.32) 

where a star denotes the complex conjugate, and the 1/2 prefactor comes from the fact that 
in Eq. (3.8) we dehned Cyy as the covariance amplitude of each independent component Wa- 
We then proceed as before, decomposing the expectation value {WeQ{v)W*^{w)) as a sum 
over all possible realizations of the SC. Since is nonzero only in holes, we fully decompose 
each realization into FL and hole elementary paths (rather that {FL+hole} sets as before). 
In the average, only the realizations such that v and w belong to the same hole H contribute 
to the net result. Hence the analogue of Eq. (6.26) is 

{\Wes{v)Wes{w)\) = {I - d) j dU p{U\v ^ Hn) Pvoh{w G Hn\v G Hn,n) \ Wes(U)\^, 

(6.33) 

where 11 is now the couple {B, R^) characterising a hole H. The (1 — d) prefactor corresponds 
to the probability that the elementary interval to which belong u is a hole. The involved 
probabilities are formally identical to the Ricci case, except that the interval length is now 
Auh instead of Av = Auh + Aufl- The integral to calculate is therefore 


{Wes{v)y^es{w)\) = (1 - a){dT^GpQiJ’f [ dll ^(11) —— 0(Auh -\v-w\ 




3 3 \ B 1 


1 2 


(6.34) 


The final result, after integration over 11 and w, is 



(6.35) 


Like in Eq. (6.31), angle brackets denote here averages with respect to the statistical properties 
of the matter clumps. 

A comparison of the covariance amplitudes and C-^, calculated with the setup and 
numerical values listed in § 6.3.1, is depicted in Fig. 2. It is clear here that Weyl covariance 
dominates over Ricci covariance. This result is characteristic of the Einstein-Straus SC model, 
where the local matter density experienced by light oscillates between p (cheese) and 0 (holes); 
this highly underestimates the fluctuations of Ricci focusing compared to reality. 
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Figure 2: Ratio C^jCyfr between the covariance amplitudes of the Ricci lensing and Weyl 
lensing in a Swiss-cheese model, as a function of the mean smoothness parameter a, for three 
different values of the redshift z = 2 (blue), z = \ (orange), and z = 0.1 (green). 


6.4 Results and comparison with ray tracing 

We now apply the general results derived in Sec. 5 with the expressions (6.31) and (6.35) for 
and C-^. After having discussed our expression of the average shear rate with respect to 
earlier works, we compare the predictions of our formalism for {D\) and var(ii>A) with the 
output of numerical ray-tracing simulations in a SC model. 


6.4.1 Shear rate and astrophysical parameter 


Introducing the expression (6.35) of in Eq. (5.30) yields the following formula for the 
average shear rate 


( |crp ) = [ dw 

Jo 


Doiw) 

Do{v) 


1 4 


(1 + ^)®, 


where we introduced a dimensionless astrophysical parameter 


(6.36) 


A = 


3 


(1 - a) 



(6.37) 


which encodes the statistical assumptions about the mass and compacity of the matter clumps. 
It also contains the main dependence with respect to the smoothness parameter a, since 
the integral of Eq. (6.36) is almost independent from it, as shown in Fig. 3. In terms of 
orders of magnitude, for a = 0, Aq ^ d/^o^ where g = GMjrl is the surface gravity of the 
central matter clumps. If they represent galaxies, then is typically of order unity, but it is 
potentially much larger for more compact objects (see table 1). 

Equation (6.36) is very similar to the ones obtained, e.g., by Gunn [62] or Kantowski [65] 
by different methods. Both get the same integral term, but their estimations of the astrophys¬ 
ical parameter differ with ours. In particular, Kantowski obtains® (Eq. (42) of Ref. [65]) 

3 /^2 — 2 \ 

= ^ \ A (6.38) 

Ho (rs) 

®Dyer and Roeder also obtained the same result, given in Eq. (25) of Ref. [69] with no derivation, but 
referring to Dyer’s PhD thesis [68]. 
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Nature of the clumps 

M 

rs 

Tc 

Ao 

galaxy clusters 

lO^^M© 

100 pc 

lOMpc 

10“^ 

galaxies 

lO^^M© 

10“^ pc 

lOkpc 

1 

stars 

Mq 

km 

10® km 

lOio 


Table 1: Typical orders of magnitude for the mass M, Schwazschild radius rg, and physical 
size Tc of three possible types of matter clumps modelled in a SC model, with the associated 
astrophysical parameter rsKH^rl) for a = 0. 


which only differs from Eq. (6.37) by the powers of rg in the averages. In a SC model where 
all the holes are identical, we thus have Ak = A, but if their masses are distributed according 
to the same distribution as in Ref. [71], then Ak/A = 1.9. Although the calculation leading 
to Eq. (6.38) is not fully detailed in Ref. [65], its discrepancy with our result (6.37) may 
be due to different statistical assumptions. In particular, we suspect that Kantowski took 
into account that bigger SC holes have a larger probability to be encountered by a light 
beam, whereas we did not—in our approach, holes are randomly placed on the line of sight, 
irrespective of their sizes. While the former is relevant in an exact SC model, the latter may 
better correspond to the actual small-scale structure of the Universe. 



redshift z 

Figure 3: Evolution of the integral of Eq. (6.36), as a function of the redshift, for three 
different smoothness parameters a = 0,0.5,1. 

6.4.2 A post-Kantowfski-Dyer-Roeder approximation 

In § 5.2.1 we have derived the general expression (5.28) of the correction = {{Da)—Dq)/Dq 
to the mean angular distance with respect to the zero-shear distance Dq —here given by the 
KDR approximation. With the formula (6.35) for in a SC model, this post-Kantowski- 
Dyer-Roeder (pKDR) term reads 


^(1) _ ^ r r n dz3 


2 

“'"io E{zi) Jo E{z 2 ) Jo E{z3) 

Do{zi)Do{z2) 



(6.39) 


with E{z) = H{z)/Ho = -yQmo(l + and where Do{z) = (1 -|- z)Dq{z) is sometimes 

called the corrected luminosity distance, here associated with the KDR distance Dq. 
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Figure 4 represents <5^^ as a function of the smoothness parameter a (4a) and of the 
redshift z (4b), comparing our calculation with the earlier result of Kantowski [65]. On Fig. 4a 
are also plotted the results of ray-tracing simulations in SC model, as described in Ref. [71]. 
Each square represents the average of (Da — Dq)/Dq over 1000 runs. These numerical results 
are thus in excellent agreement with the predictions of the stochastic leasing calculations, 
which proves its efficiency. 
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mean smoothness a 

(a) Post-KDR correction as a function of 
the smoothness parameter a, at redshift z = 1. 
Black squares are results from simulations. 


redshift z 

(b) Post-KDR correction as a function of 
redshift z for three different smoothness param¬ 
eters d = 0,0.5, 0.9. 


Figure 4: pKDR correction on the angular diameter distance = ((Da) — Dq)/Dq, at 
linear order in Weyl leasing, in SC models made of galaxylike clumps, with = 0.5. Solid 
lines correspond to our calculations and dashed lines to Kantowski’s. 


The results depicted on Fig. 4, namely ~ 10 confirm that the KDR approximation 
provides a very good effective description of the angular distance-redshift relation in SC 
models [71], at least when galaxy-like clumps are at stake. Nevertheless, since <5^^ oc A, this 
pKDR correction can become very large as the clumps are more compacts; the orders of 
magnitude given in table 1 suggests that for a SC model made of stars, 5^ ~ 10^. This 
unreasonably large number is a hint that our calculations may break down if too small 
deflectors are involved. In particular, the infinitesimal light beam approximation—on which 
both the Jacobi matrix and optical scalar formalisms are based—is not valid for describing the 
lensing of a star at cosmological distances, which rather requires a microlensing description. 
See also a discussion by Gunn in Ref. [63] on this issue. 


6.4.3 Dispersion of the angular distance 


The general equation governing the variance of the angular distance, var(IlA), based on 
second-order calculations in Weyl lensing, has been derived in § 5.3. In terms of the redshift, 
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using both Eqs. (5.14), (6.7), it reads 


6 \ d^ 

dz^ ~*~l + zy dz'^ 


H" (H'V S H' 6 

^ l + zH^ {l + zf 


1 d 2(^)' 4C^-2C^| 

(1 + zYm\ dz (1 + zYH^ (1 + zfm I 


2 Dic^ 6 r dzi r _2d^c^Y 

(1 + z)^H3 + (1 + z) 6773£,4 (1 + Zi)2i7E)2 [J^ ^2 ^ 


(6.40) 


where a prime denotes here a derivative with respect to z. To our knowledge, it is the first 
time that such a theoretical prediction of the dispersion of the angular distance though a SC 
model is proposed. This equation is solved numerically, using {^) = —(3/2)77QQ;nmo(l + 
and the expressions for and C-^ derived previously; the output is shown in Fig. 5 with, 
on Fig. 5a, a comparison with simulated data. 



mean smoothness a 

(a) Standard deviation of Da at redshift z = 1 as 
a fnnction of smoothness a. Black squares result 
from ray-tracing simulations and lines from the 
numerical integration of Eq. (6.40). 



redshift z 

(b) Standard deviation of Da as a function of 
redshift z, for three different smoothness param¬ 
eters, d = 0, 0.5,0.9. 


Figure 5: Standard deviation a^A = Y^var(T)A) of the angular distance Da in SC models, 
normalised by the KDR distance Do, as a function of the smoothness parameter a and 
redshift z. 


We see that, contrary to its average {Da), the standard deviation cda of the angular 
distance predicted by the stochastic lensing formalism does not fit with the results of ray¬ 
tracing simulations. They differ here by a factor 1.7 for d = 0. We performed a number of 
consistency checks on both the analytical and numerical sides, and found no errors. It turns 
out that such a discrepancy between theoretical and numerical results is actually a genuine 
limitation of our formalism, due to the fact that we modelled Weyl fluctuation by a Gaussian 
noise. 

Let us first show that the problem indeed comes from Weyl lensing. Formally, Fq. (6.40) 

reads 

var(i4A) = S^ + S:^, (6.41) 

where is a linear differential operator, and S^, are source functions respectively due to 
Ricci and Weyl lensing. Contrary to C^JC-p-, the ratio Sy(/ j is not necessarily small here; 


- 30 - 





















in fact, it is of order unity in the SC model used to generate the results of Fig. 5. A way to 
tune this ratio—and thus to decide which among Ricci and Weyl fluctuations dominates the 
dispersion of Da— consists in changing the lower cutoff 6niin = 'i"c of impact parameters in 
the holes. By virtue of Eq. (6.35), decreasing Tc, i.e. enhancing the compacity of the central 
clumps, increases C^. 

In Fig. 6, we compare again the predictions of the stochastic lensing formalism with 
ray-tracing results, but for two different classes of SC models: with less compact clumps 
(twice larger for the same mass, left panel); or more compact clumps (twice smaller for the 
same mass, right panel) than before. We see that the agreement between theory and numerics 
is now excellent in the first case, where S> S^, while it is slightly worse than in Fig. 5a in 
the second case, where on the contrary <C S^. The very good agreement regarding 
in both cases confirms that there are no mistakes in the evaluation of 

Such results suggest that our modelling of Ricci lensing fluctuations is more accurate 
than the one of Weyl lensing fluctuations. The weakness does not seem to be related with 
the (i-correlation hypothesis, because (i) the numerical SC model is constructed so that the 
properties of two different holes are indeed independent; (ii) the size of the holes is much 
smaller than the typical evolution scale of and (iii) this hypothesis equally applies to 
both Ricci and Weyl fluctuations, any deviation from it would therefore be manifest for any 
value of S'^/5#', which is not what we observe. 

The Gaussian hypothesis is more questionable. In the standard Langevin description of 
Brownian motion, the Gaussianity of the random force is justified by the central-limit theorem: 
during a mesoscopic time interval At, the Brownian particle is hit by many molecules, and 
the associated microscopic momentum transfers dpmicro sum into an effective transfer Ap, 
whose PDF is therefore well approximated by a Gaussian, whatever the PDF of each dpmicro- 

However, while a microscopic Brownian particle undergoes ~ 10^*^ collisions per second, a 
typical light beam in a SC models encounters only ~ 10^ holes from the source to the observer. 
The convergence towards central limit must therefore be very efficient for the Gaussian model 
to be adapted. In the case of Ricci lensing, ^ simply oscillates between 0 and the sum 
of such a random variable converges quite quickly towards the Gaussian limit, in particular 
because its support is compact. The case of Weyl lensing is different. One can easily check 
from the statistical assumptions of § 6.3.1 and the expression (6.18) of that its PDF reads 

" 3 ^ " 1 ) ' 

with Wram = 47rG/9o(l + zf' and l^ax = + 2(rh/rc)^]/3 3> l^in- This PDF thus has a 

very long algebraic tail, which drastically slows down the convergence towards central limit. 
This argument is, in our opinion, the most probable explanation of the discrepancy between 
theory and numerics observed in Fig. 5a, and of its disappearance in the left panel of Fig. 6, 
where Ricci lensing dominates. This argument shall be reinforced by the results of the next 
section. 

7 Numerical integration of the Langevin equation 

The FPK equation, Eq. (4.14) for the Jacobi matrix or Eq. (4.20) for the optical scalars, 
contains all the information necessary to characterise the statistical properties of stochastic 
lensing, provided this part of lensing can be well approximated by a Gaussian, uncorrelated 
noise (white noise). However, as mentioned before, it is in general impossible to solve explicitly 
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mean smoothness a 


mean smoothness a 




Figure 6: pKDR correction to the mean angular distance (top); normalised standard 
deviation (JD/^/Dq (middle); and ratio S^j between the Ricci and Weyl sources of variance 
(bottom), as a function of the mean smoothness parameter d at redshift z = 1, for SC models 
with two different clump densities: pc/8 2rc (left panel) and 8pc (right panel), 

where p^ = 3.5 x 10“^^ kg/m^ is the density used for the previous plots 4a, 5a. For a given SC 
hole, the minimal impact parameter h-oiin = Vc is thus respectively increased or reduced by a 
factor 2 with respect to the previous calculations. As before, squares correspond to the output 
of ray-tracing simulations, while lines are the predictions of the stochastic leasing formalism. 


the FPK equation, and one ought to rely on numerical methods to extract the statistical 
information available. From the numerical point of view, solving a partial differential equation 
is harder than tackling an ordinary differential equation and therefore, it is certainly better to 
concentrate on the Langevin equation rather than on the FPK equation. In this section, we 
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aim at solving the Langevin equation for the Jacobi matrix, Eq. (3.15) for a double purpose. 
First, we wish to show that, in the approximation of a white noise for the Ricci and Weyl 
lensing, the ray-tracing and analytical results of the previous section are well re-produced by 
directly solving the Langevin equation. Second, we would like to probe the effect of relaxing 
the Gaussian approximation: we will show that in the SC model, if the ‘true’ PDF of 
Eq. (6.42), is used, the discrepancy in the fluctuations of Da between the ray-tracing results 
and the analytical estimates coming from the FPK equation can clearly be attributed to the 
non-Gaussianity of the noise and the lack of convergence towards the central limit when the 
number of holes encountered is too small. 

7.1 The stochastic Euler method 

We begin by a short exposition of the numerical discretisation of the general Langevin 
equation (4.1). For an infinitesimal time step df, we can rewrite this equation as 

X{t + dt) = X{t) + f{X,t) dt + L{X, t)N{t) dt. (7.1) 

Noting dB{t) = N{t)dt, this becomes simply 

X{t + dt) = X{t) + /(X, t)dt + L{X, t) dB{t), (7.2) 

which, after discretisation, gives the Euler approximation to the Langevin equation 

X{ti+i) = X{ti) + f[X{ti),ti] At + L[X{ti),ti] AB{ti), (7.3) 

where we have assumed, for simplicity, a constant time step At. As discussed in Sec. 4.1, 
if the noise At is a white noise, then S is a Brownian motion, i.e. its increment AB is a 
zero-mean Gaussian process with covariance matrix 

(AB{ti)AB'^{tj)) = Q{ti)6ij At (no summation over i). (7.4) 

In practice, simulating one realisation of the process X (t) is thus identical to numerically 
solving an ordinary differential equation, except that at each time step tj the quantity AB{ti) 
is randomly picked, according to a Gaussian PDF with variance Q(ti)At, and independently of 
the other steps. This means that the components of the stochastic term AB have fluctuations 
on the order of \/At: the stochastic Euler method only converges as \/At, instead of At for 
its deterministic counterpart. This limitation will not be a problem in what follows. 

Note also that we can still apply this discretisation if the noise is not Gaussian, but the 
term N(t)At must then be evaluated from the true PDF of N. The main caveat, in this case, 
lies on the fact that the simulation is no longer resolution independent (see § 7.3). 

7.2 Application to the Swiss-cheese model — Gaussian case 

Let us now turn to our specific Langevin equation for the Jacobi matrix (4.14), in the SC 
model. Using the relationship between affine parameter and redshift, we can rewrite it as 

^ = ( 7 . 5 ) 

Discretising this equation with a constant redshift step Az (for simplicity), one gets 

Jfc+i = Jk + [M{zk)JkAz + Ljac(Tfc)AS(2:fc)] , (7.6) 
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with Jk = J{zk), and where the covariance matrix of AB involves the diffusion matrix of 
Eq. (3.10) according to 

{AB{zk)AB^{zk)) = Q{zk)Az (7.7) 

= diag(C'^, C^, C^)Az. (7.8) 

At each time step, the quantity AB^ = [AB^, AB^^, AB^^) is obtained by randomly picking 
ABgi and AB^^ according to a zero-mean Gaussian distribution with variance C^Az and 
respectively. 

Using the expressions for C^{z) and C-^{z) found in the case of the SC model, Eqs. (6.31) 
and (6.35) respectively, we can now integrate numerically the Langevin equation. Results 
are presented in Eig. 7, for the same set of parameters as those used in the previous section, 
i.e. with a standard distribution of galaxy-type holes. Statistical averages are performed 
over 1000 realisations of J{z), for each possible value of d, each realisation being simulated 
according to the stochastic Euler method with a redshift step Az = 10“^. The agreement 
with the results from the FPK approach is striking, and provides strong support for the 
analytical expressions found previously. Compared to ray-tracing simulations, the pKDR 
corrections to (Ha) are very accurately reproduced, but the dispersion of Ha suffers from the 
same systematic underestimation as in the previous section. If one could resolve this tension, 
because the numerical integration of the Langevin equation is much faster than ray-tracing 
simulations, it would provide an efficient way to estimate statistical quantities. 




Figure 7: Results of the numerical integration of the Langevin equation with a Gaussian 
noise (empty purple circles), compared with analytical calculations (blue lines) and ray-tracing 
simulations (black squares) in the SC model. Left panel (same as Fig. 4a): pKDR correction, 
in percent, to the angular distance at z = 1, as a function of a. Right panel (same as Fig. 5a): 
fractional dispersion, in percent, of Ha at z = 1 as a function of a. 


7.3 Beyond the Gaussian approximation 

Now that we have shown that numerically integrating the Langevin equation leads to the same 
results as the use of the FPK equation in the Gaussian noise limit, we would like to show that 
the discrepancy between these results and the ray-tracing results stems from non-Gaussianity. 
Indeed, as discussed in the previous section, the Weyl lensing is very poorly described by 
a Gaussian noise, since its actual PDF in a SC model presents a long non-Gaussian tail, 
corresponding to not so rare events during which light rays pass very close to masses and 
experience significant tidal distortions. In order to probe this effect, in this subsection, we 
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limit ourselves to the case a = 0 in which the Ricci lensing is zero, thus isolating effects due 
to a pure Weyl lensing. We also use an SC model with one size of holes for simplicity. 

We come back to Eq. (7.5) but we no longer treat the noise term as the increments AB 
of a Brownian motion B, and replace it by AB such that 

AB^^ = l^effl sin 2/3, (7.10) 

where fi is uniformly distributed within [0,27r], while the PDF of |#^fr| is given by Eq. (6.42). 
One can check that the above choice ensures that (^ABAB"^) = d\a.g{C^,Cyy,C^)Az as 
before. In other words, the resulting AB is a non-Gaussian process whose first two moments 
match the ones of the Gaussian model. Note that this is somehow artificial, because the 
noise modelled by AB now depends on the resolution Az used for integrating the Langevin 
equation. This can be understood as follows. Suppose one solves the Langevin equation 
with two different resolutions: a low resolution (LR) A^lR; and a high resolution (HR) 
Azhr = Azlr/u. During a given interval [z, z + Azlr], the HR simulation performs n steps, 
and the effective noise associated with the set of these n steps reads 


n— 1 

A.Bhr(2; —)■ 2 ; + Azlr) = ^ ABnKizk Zk+i), (7-11) 

fc =0 


with Zk = z + /cAzhr- Contrary to the Gaussian case, the above sum is not equal to ABlR) 
because any random variable does not enjoy the invariance under addition; in particular, for 
n —>• 00 it becomes Gaussian itself, by virtue of the central limit theorem. 

We therefore expect the output of numerical integration of the Langevin equation with 
AS (i) to depend on the resolution A 2 ;, and (ii) to converge towards the Gaussian case for 
A 2 ; —)• 0. This is illustrated in Fig. 8, where the mean and dispersion of the angular distance 
at z = 1, obtained by integrating Langevin equation in the Gaussian and non-Gaussian 
cases, are plotted as a function of Az. We also indicate, for comparison, the analytical and 
ray-tracing results. A number of comments shall be formulated about those figures. First, all 
the results on the mean angular distance (Da)— more precisely, its pKDR correction —are 
in excellent agreement. It is not the case concerning the dispersion crop, of Da- Then, as 
expected, the Gaussian numerical results coincide with the analytical calculations, as well 
as the non-Gaussian result for Az —)• 0. The latter however depart from the formers as Az 
increases, and coincides with the ray-tracing results for Az ~ 2.5 x 10“^. This particular value 
can be understood as follows: physically speaking, a non-Gaussian Langevin simulation with 
redshift step Az corresponds to a SG model where successive holes are typically separated 
by Az, that is z/N where z is the redshift of the source and N the typical number of holes 
between the source and the observer. As a matter of fact, with the parameters used for 
generating Fig. 8, the average number of holes encountered by a photon is on the order of 
3000, corresponding to a Az ~ 3 x 10““^, which is very close to the value 2.5 x 10“^ where 
the ray-tracing and the stochastic non-Gaussian results match. 

This confirms our point that, in the SG models investigated here, the typical number 
of collisions is marginally too small to warrant a treatment of the lensing in terms of a pure 
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white noise, i.e. with a FPK equation. This understanding of the problem provides two ways 
of escaping from it: (1) dealing with smaller-scale structures; (2) increasing the redshift z of 
the source. In both situation, the number N of deflectors, that is the physical resolution of 
the problem, is increased, which should thus improve the agreement between exact ray-tracing 
results and the analytical FPK calculations. A quantitative criterion, allowing us to estimate 
the precision of the FPK approach, remains nevertheless to be determined. 
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Figure 8: pKDR correction to {D\) (top panel) and dispersion of Da (bottom panel) at 
z = 1, computed from numerical integration of the Langevin equation (7.5) with a Gaussian 
noise AS (empty circles) or a non-Gaussian noise AS (filled circles), as functions of the 
redshift step Az used in the Euler method. Each circle is obtained from the statistical 
properties of a sample of 1000 realisations of J. Eor comparison, dotted lines indicate the 
output of ray-tracing simulations, while the solid lines are the analytical predictions of the 
stochastic leasing formalism, i.e. given by Eqs. (6.39), (6.40). The smoothness parameter of 
the underlying SG model is d = 0; all holes have the same mass M = IO^^Mq and density pc- 
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8 Conclusion 


In this article, we proposed a new theoretical framework in which the gravitational leasing 
caused by the small-scale structure of the Universe is treated as a diffusion process. The Sachs 
equations governing the propagation of narrow light beams were provided with stochastic 
components modelled as white noises. We derived the associated Fokker-Planck-Kolmogorov 
equations for the PDF of the Jacobi matrix and the optical scalars. We used them to deduce 
(1) the corrections to the mean angular distance due to Weyl lensing, and (2) a differential 
equation for the dispersion of the angular distance. These results depend on three free 
functions, namely the mean Ricci lensing (^) and the covariance amplitudes of Ricci and 
Weyl lensing C^,C#', which need to be specified from a model. As both an illustration 
and a test of this new formalism, we applied it to Einstein-Straus Swiss-cheese models. The 
results on (Da) offer an extension to the Kantowski-Dyer-Roeder approximation, in excellent 
agreement with numerical simulations. The theoretical predictions for the variance of Da 
are however systematically lower than their numerical counterpart. We located the origin of 
this discrepancy in the actual non-Gaussianity of Weyl lensing, which cannot be captured by 
the FPK approach. This was confirmed by direct simulations of the Langevin equation with 
Gaussian and non-Gaussian source terms. 

This new approach has the advantage of dealing with small-scale lensing in a mathemat¬ 
ically consistent and efficient way, without the need for computationally expensive ray-tracing 
simulations. It complements the standard description of weak lensing caused by the large-scale 
structure, allowing for the effect of smaller scales. It is also very flexible, in the sense that it 
can be applied, in principle, to any model for the distribution of matter on those scales. 

The main limitation of our formalism, under its present form, lies in the assumption of 
Gaussianity. This hypothesis is indeed central in the general derivation of the Fokker-Planck- 
Kolmogorov equation, on which our main results are based, but it may not hold in the actual 
Universe, as illustrated on the particular example of Swiss-cheese models. There is, on the 
mathematics side, active ongoing research on stochastic processes with non-Gaussian noises. 
Unfortunately no dehnite standard prescription about how to modify the FPK equation has 
been established so far, which is the reason why we did not enter into such discussions in the 
present article. Nevertheless, from a practical point of view, we empirically checked that the 
Gaussian limit provides good estimations of the lensing quantities, as far as only their mean 
and variance are concerned: the largest discrepancies are expected to appear for higher-order 
moments. 

In the future, we plan to apply the stochastic lensing framework to more realistic models 
than the Swiss cheese, in particular for comparing its output with the numerical results of 
Refs. [86, 87]. We also intend to include the effect of peculiar velocities, which is not expected 
to contain major difficulties. On longer terms, we aim at explicitly combining our formalism 
with the standard perturbation theory, which would ideally provide a consistent multiscale 
treatment of cosmological lensing. This will however require to establish a quantitative 
criterion about the transition scale from one behaviour to the other. Finally, we emphasize 
the very general character of the approach presented here, which may also be applied to 
spacetimes with very different symmetries, to treat e.g. the microlensing due to stars in a 
galaxy, or the effect of a stochastic background of gravitational waves. 
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A Geometric optics in curved spacetime 

This appendix summarises textbook elements about the propagation of light in arbitrary 
spacetimes, which aims at supplementing the relatively sharp presentation of Sec. 2. Eor more 
general introductions, see Refs. [92, 111-114]. 

A.l Description of a light beam 

A light beam is a collection of light rays, that is, a bundle of null geodesics {v i—>■ x^{v,y°')} 
converging at a given event (here taken to be the observation event O), where the two 
coordinates (y“)a=i ,2 label the rays, while v is the affine parameter along them. There is no 
need for a fourth coordinate because the beam entirely belongs to the lightcone of O, which is 
an isophase hypersurface. 

The wave four-vector = dx^/dv is a null vector field, tangent to the rays = cst. It 

satisfies the null geodesic equations 

= 0, and k^Vuk^ = 0. (A.l) 

Besides, the relative behaviour of two neighbouring geodesics of the bundle, x^(-,y“) 
and a:^(-,y“ -|- hy“), is described by their connecting vector = [dx^/dy°‘)5y°'. If the origin 
u = 0 of the affine parametrisation of all rays is taken at O, then 

k>^^^. = 0. (A.2) 

As soon as the condition (A.2) is satisfied, the evolution of along the beam is governed by 
the geodesic deviation equation 

k^kf^VaVfse = (A.3) 

where R^uap is the Riemann tensor. 

A.2 The Sachs formalism 

Consider an observer, with four-velocity (u^u^ = —1), who crosses the light beam. The 
spatial direction of propagation of the beam, relative to this observer, is defined as the opposite 
of the direction in which the observer must look to detect a signal. It is spanned by a purely 
spatial unit vector 

= 0, = 1, (A.4) 

such that (remember that we took k^ future oriented in this article) 

F =-w(u^ + d^), (A.5) 
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where 


uj = 27rz^ = Ufj,k^ (A.6) 

is the cyclic frequency of the light signal in the observer’s rest frame. Note that d£ = Lodv is 
the proper distance (measured by the observer) travelled by light for a change du of the affine 
parameter. The redshift 2 ; is defined as the relative change between the emitted frequency Vg,, 
in the source’s frame, and the observed frequency in the observer’s frame, that is 


14 _ u^kf,{vg) 
U^kfj,{Vo) ' 


(A.7) 


Suppose that the observer measures the size and the shape of the light beam. For that 
purpose, he must use (and thus define) a (spatial) screen orthogonal to the line of sight. This 
screen is spanned by the so-called Sachs basis (s(^)Ae{i, 2 }! defined by 


— ^ABi 


(A.8) 


and by the transport property (A.9) below. The projections indicate the relative 

position, on the observer’s screen, of the light points corresponding to two neighbouring rays 
separated by Thus, it encodes all the information about the size and shape of the beam. 

Consider a family of observers u^{v), along the beam, who wants to follow the evolution 
of the shape of the beam (typically for shear measurements). For that purpose, they must all 
use the “same” Sachs basis, in order to avoid any spurious rotation of the pattern observed on 
the screens. This is ensured by imposing that the Sachs basis is a parallel transported as 


S^,ukPVps\ = 0, (A.9) 

where 

= 5^^s^Xs'^B = g^''+ -dP(f (A.IO) 

is the screen projector. The reason why cannot be completely parallel-transported is that, 
in general, is not'. 

The evolution of ^a, with light propagation, is determined by projecting the geodesic 
deviation equation (A.3) on the Sachs basis. The result is known as the Sachs equation [112, 
115, 116], 

(A.ll) 

where 

= Rp.a/sk’^k^s'Xs^B (A-12) 

is the screen-projected Riemann tensor, usually called the optical tidal matrix. The properties 
of the Riemann tensor imply that this matrix is symmetric, TZab = R-ba- Note that the 
position of the screen indices (A, B,...) does not matter, since they are raised and lowered 
by Sab- In tbis article, to alleviate the notation, we use bold symbols for quantities with 
screen indices, and an overdot for derivatives with respect to the affine parameter v. The 
Sachs equation (A.ll) thus becomes 

i = (A.13) 

^In fact, it is also possible to choose a family of observers such that the four-velocity field is parallel- 

transported along the beam, without affecting the optical equations [112]. In this case, however, the observers 

are generally not comoving, and thus have no clear cosmological interpretation. 
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A.3 Evolution in terms of potentials 

It is interesting to note that the Sachs equation can be reformulated in terms of a potential, as 


Ia = 


dV 


(A.14) 


with 

V{T>) = -^^^Uab^b = (A.15) 

This equally applies to the Jacobi matrix equation Vab = —dVji^cld'DAB, with Vjac(2^) = 
-{l/2)VABnAc'DcB = -(i/2)tr (p'^nv). 

Regarding the optical scalars, the Sachs equations can be rewritten as® 


dv 



(g^)Rscal + 



(A.17) 


where 


f^cal — 


03 

3 




(A.18) 


and the focusing scalars and W are here treated as an external force; they could also have 
been included in the potential according to 


t4cal = Rscal-^0-^fT. (A. 19) 

Note however that I4cai does not depend explicitly on u, while Rgcai and Vjac generally do, 
because of the presence of ^ and W. 


A.4 Decompositions of the Jacobi matrix 

As defined in Sec. 2, the Jacobi matrix relates the physical separation ^^{v) of two neighbouring 
rays of a beam at v to their observed separation ^^(0), as 

i^{v) = V^s{v)i^{^)- (A.20) 

When applied to an observed image, 'D{vs) thus returns the intrinsic physical properties of 
the source. 


A.4.1 General decomposition 

As any 2x2 (nonsymmetric) matrix, the Jacobi matrix has 4 degrees of freedom. It can be 
decomposed in a way that highlights the geometrical transformations between the source and 
the image. First of all, up to frequency factor uJq fixed to 1 in this article, the determinant of 
the Jacobi matrix is related to the angular diameter distance as 


Dliv) 


area of the source at u 

observed angular size 


(A.21) 


®As usual, we define the complex derivative as 


df 1 f df .df\ 

da~2\dai 


(A.16) 


for a — ai + ia 2 . 
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Factorising the determinant, we are left with a 2 x 2 matrix of determinant 1, which can be 
decomposed as the product between a symmetric matrix and the exponential of a symmetric 
traceless matrix: 

■D = dJ '''') ■ (A.22) 

\—simpcosipJ V 72 7 i/ 

The exponential matrix can also be diagonalised by defining 7 > 0 and as 

( 7 I) 72 ) = 7 (cos 2 i?, - sin 2 i?), (A.23) 


so that 


V = Da 


cos 'll: — sin 'll: 
sincos 'll: 


cos 1 ? — sin /e 0 
sin?? cos?? / V 0 e"^ 


cos '& sin ?? 
— sin ?? cos ?? 


(A.24) 


This decomposition shows that, in order to reconstruct the physical properties of a light 
source from its observed image, one must: 


1 . Contract it by a factor e~'^ along a direction inclined of ?? with respect to the Sachs 
basis, and stretch it by a factor e"^ along the orthogonal direction. This represent the 
net shear, which preserves the area of the image. 

2. Rotate anticlockwise the result by an angle 'll:. 

3. Scale it with Da to turn angles into lengths. 

Note that, by virtue of Etherington’s reciprocity relation [92, 117], which stipulates that 
the Jacobi matrix obtained by integrating the Sachs equation (A. 11 ) from the observer O to 
the source S or from the source to the observer are opposite and transposed with respect to 
each other,® 

T>{S ^0) = ^ S), (A.26) 

the net shear 7 is independent of the sense in which this integration is made. This is the general 
nonperturbative generalization of the shear reciprocity relation mentioned in appendix A of 
Ref. [39]. 


A.4.2 Perturbative case 

Usually, when dealing with weak leasing as caused by perturbations with respect to Minkowski 
or Friedmann-Lemaitre spacetimes, one uses that at background level both shear and rotation 
vanish so that D = Da 1- The decomposition (A.24) can then be expanded at first order in 
7 i, 72 and ?/; to get the definition of the amplifieation matrix as 

T> = At> + 0(2) (A.27) 


with 


A-fi:-7i 72+ V’ ^ 

\ 72- V’ I- K + Ji) 


(A.28) 


® This can directly be shown from the fact that = 7 ?,, which implies that, for any two vi,V 2 the 
function C(v) defined by C(v) = 'D (v <— vi)T>(v <— V 2 ) - 'D'^iv vi)'D{v •<— V 2 ) is a constant. Writing 
C(vi) = C{v 2 ), we conclude that 

P(ni ^ V 2 ) = -'D'^{v 2 ^ vi). (A.25) 

This relation is central to the derivation of the distance duality relation Dl = {\ + z)^D a. The latter is however 
more easily achieved by abandoning the convention oiq = 1. See e.g. Ref. [118] for further details. 
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where the convergence is defined as 

K = l- hTA= +0(2). (A.29) 

2 Da 

The rotation angle "0 can be proved to be on the order of 7 ^, and can thus be omitted at linear 
order, which yields the standard form of the amplification matrix. The decomposition (A.24) 
is however much more relevant in nonperturbative cases, or when either the shear or the 
rotation does not vanish at background level; see e.g. Ref [119]. 

References 

[1] P. Peter and J.-P. Uzan, Primordial Cosmology. Oxford University Press, 2009. 

[2] S. R. Green and R. M. Wald, How well is our Universe described by an FLRW model?, Classical 
and Quantum Cravity 31 (Dec., 2014) 234003, [arXiv; 1407.8084]. 

[3] T. Buchert, M. Carfora, G. F. R. Ellis, E. W. Kolb, M. A. H. MacCallum, J. J. Ostrowski, 

S. Rasanen, B. F. Roukema, L. Andersson, A. A. Goley, and D. L. Wiltshire, Is there proof that 
backreaction of inhomogeneities is irrelevant in cosmology?, ArXiv e-prints (May, 2015) 

[arXiv: 1505.07800]. 

[4] S. R. Green and R. M. Wald, Comments on Backreaction, ArXiv e-prints (June, 2015) 

[arXiv: 1506.06452]. 

[5] T. Futamase and M. Sasaki, Light Propagation and the Distance Redshift Relation in a Realistic 
Inhomogeneous Universe, Phys. Rev. D40 (1989) 2502. 

[6] A. Cooray, D. Holz, and D. Huterer, Cosmology from supernova magnification maps, Astrophys. 
J. 637 (2006) L77-L80, [astro-ph/0509579]. 

[7] S. Dodelson and A. Vallinotto, Learning from the scatter in type ia supernovae, Phys. Rev. D74 
(2006) 063515, [astro-ph/0511086]. 

[8] P. Valageas, Weak gravitational tensing effects on the determination of omega-0 and lambda 
from sneia, Astron. Astrophys. 354 (2000) 767, [astro-ph/9904300]. 

[9] G. Bonvin, R. Durrer, and M. A. Gasparini, Fluctuations of the luminosity distance, Phys. Rev. 
D73 (2006) 023523, [astro-ph./0511183]. [Erratum: Phys. Rev.D85,029901(2012)]. 

[10] N. Meures and M. Bruni, Redshift and distances in a Lambda-CDM cosmology with non-linear 
inhomogeneities, Mon. Not. Roy. Astron. Soc. 419 (2012) 1937, [arXiv: 1107.4433]. 

[11] I. Ben-Dayan, M. Gasperini, G. Marozzi, E. Nugier, and G. Veneziano, Average and dispersion 
of the luminosity-redshift relation in the concordance model, JCAP 6 (June, 2013) 2, 

[arXiv: 1302. 0740]. 

[12] I. Ben-Dayan, M. Gasperini, G. Marozzi, E. Nugier, and G. Veneziano, Do Stochastic 
Inhomogeneities Affect Dark-Energy Precision Measurements?, Physical Review Letters 110 
(Jan., 2013) 021301, [arXiv : 1207. 1286]. 

[13] O. Umeh, G. Clarkson, and R. Maartens, Nonlinear relativistic corrections to cosmological 
distances, redshift and gravitational leasing magnification: 1. Key results. Classical and 
Quantum Gravity 31 (Oct., 2014) 202001, [arXiv : 1207.2109]. 

[14] O. Umeh, C. Clarkson, and R. Maartens, Nonlinear relativistic corrections to cosmological 
distances, redshift and gravitational leasing magnification: 11. Derivation, Classical and 
Quantum Gravity 31 (Oct., 2014) 205001, [arXiv : 1402.1933]. 

[15] E. Barausse, S. Matarrese, and A. Riotto, Effect of inhomogeneities on the luminosity 
distance-redshift relation: Is dark energy necessary in a perturbed universe ?, Phys. Rev. D 71 
(Mar., 2005) 063537, [astro-ph/0501152]. 


- 42 - 



[16] I. Ben-Dayan, G. Marozzi, F. Nugier, and G. Veneziano, The second-order luminosity-redshift 
relation in a generic inhomogeneous cosmology, JCAP 11 (Nov., 2012) 45, [arXiv; 1209.4326]. 

[17] S. Andrianomena, G. Clarkson, P. Patel, O. Umeh, and J.-P. Uzan, Non-linear relativistic 
contributions to the cosmological weak-lensing convergence, JCAP 1406 (2014) 023, 

[arXiv: 1402.4350]. 

[18] C. Clarkson, O. Umeh, R. Maartens, and R. Durrer, What is the distance to the CMB?, JCAP 
11 (Nov., 2014) 36, [arXiv: 1405.7860]. 

[19] T. W. B. Kibble and R. Lieu, Average Magnification Effect of Clumping of Matter, ApJ 632 
(Oct., 2005) 718-726, [astro-ph/0412275]. 

[20] N. Kaiser and J. A. Peacock, On the Bias of the Distance-Redshift Relation from Gravitational 
Tensing, ArXiv e-prints (Mar., 2015) [arXiv: 1503.08506]. 

[21] C. Bonvin, C. Clarkson, R. Durrer, R. Maartens, and O. Umeh, Cosmological ensemble and 
directional averages of observables, JCAP 7 (July, 2015) 40, [arXiv: 1504.01676]. 

[22] C. Bonvin, C. Clarkson, R. Durrer, R. Maartens, and O. Umeh, Do we care about the distance 
to the CMB? Clarifying the impact of second-order lensing, JCAP Q (June, 2015) 50, 

[arXiv: 1503.07831]. 

[23] A. Einstein and E. Straus, The influence of the expansion of space on the gravitation fields 
surrounding the individual stars. Rev. Mod. Phys. 17 (1945) 120. 

[24] A. Einstein and E. Straus, Corrections and additional remarks to our paper: The influence of 
the expansion of space on the gravitation fields surrounding the individual stars. Rev. Mod. 
Phys. 18 (1945) 148. 

[25] V. Marra, E. W. Kolb, S. Matarrese, and A. Riotto, On cosmological observables in a 
swiss-cheese universe, Phys. Rev. D76 (2007) 123004, [arXiv: 0708.3622]. 

[26] N. Brouzakis, N. Tetradis, and E. Tzavara, Light Propagation and Large-Scale Inhomogeneities, 
JCAP 0804 (2008) 008, [astro-ph/0703586]. 

[27] N. Brouzakis, N. Tetradis, and E. Tzavara, The effect of large scale inhomogeneities on the 
luminosity distance, JCAP 2 (Eeb., 2007) 13, [astro-ph/0612179]. 

[28] T. Biswas and A. Notari, Swiss-Cheese Inhomogeneous Cosmology and the Dark Energy 
Problem, JCAP 0806 (2008) 021, [astro-ph/0702555]. 

[29] V. Marra, E. W. Kolb, and S. Matarrese, Light-eone averages in a Swiss-cheese universe, Phys. 
Rev. D 77 (Jan., 2008) 023003, [arXiv:0710.5505]. 

[30] T. Clifton and J. Zuntz, Hubble Diagram Dispersion Erom Large-Scale Structure, Mon. Not. 
Roy. Astron. Soc. 400 (2009) 2185, [arXiv: 0902.0726]. 

[31] S. J. Szybka, On light propagation in Swiss-Cheese cosmologies, Phys. Rev. D84 (2011) 044011, 
[arXiv: 1012. 5239]. 

[32] R. A. Vanderveld, E. E. Flanagan, and 1. Wasserman, Luminosity distance in ’Swiss cheese’ 
cosmology with randomized voids: 1. Single void size, Phys. Rev. D78 (2008) 083511, 

[arXiv: 0808. 1080]. 

[33] W. Valkenburg, Swiss Cheese and a Cheesy CMB, JCAP 0006 (2009) 010, [arXiv: 0902. 4698]. 

[34] K. Bolejko, The effect of inhomogeneities on the distance to the last scattering surface and the 
accuracy of the CMB analysis, JCAP 2 (Feb., 2011) 25, [arXiv: 1101.3338]. 

[35] E. E. Flanagan, N. Kumar, 1. Wasserman, and R. A. Vanderveld, Luminosity distance in Swiss 
cheese cosmology with randomized voids. 11. Magnification probability distributions, Phys. Rev. 
D85 (2012) 023510, [arXiv : 1109. 1873]. 


- 43 - 


[36] E. E. Flanagan, N. Kumar, and I. Wasserman, Luminosity distance in Swiss cheese cosmology 
with randomized voids and galaxy halos, Phys. Rev. D88 (2013), no. 4 043004, 

[arXiv: 1207. 3711], 

[37] K. Bolejko, C. Clarkson, R. Maartens, D. Bacon, N. Meures, and E. Beynon, Antilensing: The 
Bright Side of Voids, Physical Review Letters 110 (Jan., 2013) 021302, [arXiv; 1209.3142]. 

[38] M. Lavinto, S. Rasanen, and S. J. Szybka, Average expansion rate and light propagation in a 
cosmological Tardis spacetime, JCAP 12 (Dec., 2013) 51, [arXiv; 1308.6731]. 

[39] M. Lavinto and S. Rasanen, CMB seen through random Swiss Cheese, ArXiv e-prints (July, 
2015) [arXiv; 1507.06590]. 

[40] K. Bolejko, The Szekeres Swiss Cheese model and the CMB observations, Ceneral Relativity and 
Gravitation 41 (Aug., 2009) 1737-1755, [arXiv; 0804.1846]. 

[41] K. Bolejko and M.-N. Celerier, Szekeres Swiss-cheese model and supernova observations, Phys. 
Rev. D82 (Nov., 2010) 103510, [arXiv; 1005.2584]. 

[42] A. Peel, M. A. Troxel, and M. Ishak, Effect of inhomogeneities on high precision measurements 
of cosmological distances, Phys. Rev. D QO (Dec., 2014) 123536, [arXiv; 1408.4390]. 

[43] M. A. Troxel, M. Ishak, and A. Peel, The effects of structure anisotropy on lensing observables 
in an exact general relativistic setting for precision cosmology, JCAP 3 (Mar., 2014) 40, 
[arXiv; 131 1.5936]. 

[44] J. Adamek, E. Di Dio, R. Durrer, and M. Kunz, Distance-redshift relation in plane symmetric 
universes, Phys. Rev. D89 (2014), no. 6 063543, [arXiv; 1401. 3634]. 

[45] K. Bolejko and P. G. Ferreira, Ricci focusing, shearing, and the expansion rate in an almost 
homogeneous Universe, JCAP 1205 (2012) 003, [arXiv; 1204.0909]. 

[46] E. V. Linder, Light propagation in generalized Friedmann universes, A&A 206 (Nov., 1988) 
190-198. 

[47] E. V. Linder, Transition from Clumpy to Smooth Angular Diameter Distances, ApJ 497 (Apr., 
1998) 28-31, [astro-ph/9707349]. 

[48] E. V. Linder, Averaging Lnhomogeneous Universes: Volume, Angle, Line of Sight, ArXiv 
Astrophysics e-prints (Jan., 1998) [astro-ph/9801122]. 

[49] S. Rasanen, Light propagation in statistically homogeneous and isotropic dust universes, JCAP 
0902 (2009) oil, [arXiv;0812.2872]. 

[50] S. Rasanen, Light propagation in statistically homogeneous and isotropic universes with general 
matter content, JCAP 1003 (2010) 018, [arXiv;0912.3370]. 

[51] S. Rasanen, Light propagation and the average expansion rate in near-FRW universes, Phys. 
Rev. D85 (2012) 083528, [arXiv; 1107.1176]. 

[52] M. Gasperini, G. Marozzi, F. Nugier, and G. Veneziano, Light-cone averaging in cosmology: 
formalism and applications, JCAP 7 (July, 2011) 8, [arXiv; 1104. 1167]. 

[53] 1. Ben-Dayan, M. Gasperini, G. Marozzi, F. Nugier, and G. Veneziano, Backreaction on the 
luminosity-redshift relation from gauge invariant light-cone averaging, JCAP 4 (Apr., 2012) 36, 
[arXiv; 1202. 1247]. 

[54] S. Bagheri and D. J. Schwarz, Light propagation in the averaged universe, JCAP 10 (Oct., 
2014) 73, [arXiv; 1404.2185]. 

[55] C. Clarkson, G. F. R. Ellis, A. Faltenbacher, R. Maartens, O. Umeh, and J.-P. Uzan, 

(Mis-)Interpreting supernovae observations in a lumpy universe, Mon. Not. Roy. Astron. Soc. 
426 (2012) 1121-1136, [arXiv; 1109.2484]. 


- 44 - 


[56] Y. Zel’dovich, Observations in a Universe Homogeneous in the Mean, Sov. Astron. Lett. 8 
(1964) 13. 

[57] R. P. Feynman, 1964. Unpublished colloquium given at the California Institute of Technology. 

[58] V. M. Dashevskii and Y. B. Zel’dovich, Propagation of light in a nonhomogeneous nonflat 
universe ii, Sov. Astronom. 8 (1965) 854. 

[59] V. M. Dashevskii and V. I. Slysh, On the Propagation of Light in a Nonhomogeneous Universe, 
Astron. Zh. 42 (1965) 863. 

[60] V. M. Dashevskii and V. I. Slysh, On the Propagation of Light in a Nonhomogeneous Universe, 
Sov. Astronom. 9 (Feb., 1966) 671. 

[61] B. Bertotti, The Luminosity of Distant Galaxies, Royal Society of London Proceedings Series A 
294 (Sept., 1966) 195-207. 

[62] J. E. Gunn, On the Propagation of Light in Inhomogeneous Cosmologies. I. Mean Effects, ApJ 
150 (Dec., 1967) 737. 

[63] J. E. Gunn, A Fundamental Limitation on the Accuracy of Angular Measurements in 
Observational Cosmology, ApJ 147 (Jan., 1967) 61. 

[64] S. Weinberg, Apparent luminosities in a locally inhomogeneous universe, ApJL 208 (Aug., 1976) 
L1-L3. 

[65] R. Kantowski, Corrections in the Luminosity-Redshift Relations of the Homogeneous Friedmann 
Models, ApJ 155 (Jan., 1969) 89. 

[66] C. Dyer and R. Roeder, The Distance-Redshift Relation for Universes with no Intergalactic 
Medium, Astrophys. J. 174 (1972) L115. 

[67] C. C. Dyer and R. C. Roeder, Distance-Redshift Relations for Universes with Some Intergalactic 
Medium, ApJL 180 (Eeb., 1973) L31. 

[68] C. C. Dyer, Observational Aspects of Locally Inhomogeneous Cosmological Models. PhD thesis, 
University of Toronto (Canada), 1973. 

[69] C. C. Dyer and R. C. Roeder, Observations in Locally Inhomogeneous Cosmological Models, 
ApJ 189 (Apr., 1974) 167-176. 

[70] R. C. Roeder, Apparent magnitudes, redshifts, and inhomogeneities in the universe, ApJ 196 
(Mar., 1975) 671-673. 

[71] P. Fleury, Swiss-cheese models and the Dyer-Roeder approximation, JCAP 1406 (2014) 054, 
[arXiv: 1402. 3123]. 

[72] D. E. Holz and R. M. Wald, A New method for determining cumulative gravitational lensing 
effects in inhomogeneous universes, Phys. Rev. D58 (1998) 063501, [astro-ph/9708036]. 

[73] T. Okamura and T. Futamase, Distance-Redshift Relation in a Realistic Inhomogeneous 
Universe, Prog. Theor. Phys. 122 (2009) 511-520, [arXiv : 0905.1160]. 

[74] J.-P. Bruneton and J. Larena, Observables in a lattice Universe, Class. Quant. Grav. 30 (2013) 
025002, [arXiv: 1208.1411]. 

[75] T. Clifton and P. G. Ferreira, Archipelagian Cosmology: Dynamics and Observables in a 
Universe with Discretized Matter Content, Phys. Rev. D80 (2009) 103503, [arXiv:0907.4109]. 
[Erratum: Phys. Rev.D84,109902(2011)]. 

[76] T. Clifton and P. G. Eerreira, Errors in estimating D? due to the fluid approximation, JCAP 10 
(Oct., 2009) 26, [arXiv: 0908.4488]. 

[77] T. Clifton, P. G. Ferreira, and K. O’Donnell, Improved treatment of optics in the 
Lindquist-Wheeler models, Phys. Rev. D 85 (Jan., 2012) 023502, [arXiv : 1110.3191]. 


-45 - 


[78] R. Kantowski, The Effects of Inhomogeneities on Evaluating the Mass Parameter Vim and the 
Cosmological Constant ApJ 507 (Nov., 1998) 483-496, [astro-ph/9802208]. 

[79] A. Cooray, D. Huterer, and D. Holz, Problems with pencils: lensing covariance of supernova 
distance measurements, Phys. Rev. Lett. 96 (2006) 021301, [astro-ph/0509581]. 

[80] D. Sarkar, A. Amblard, D. E. Holz, and A. Cooray, Lensing and Supernovae: Quantifying The 
Bias on the Dark Energy Equation of State, Astrophys. J. 678 (2008) 1, [arXiv : 0710. 4143]. 

[81] P. Fleury, H. Dupuy, and J.-P. Uzan, Can all cosmological observations be accurately interpreted 
with a unique geometry?, Phys.Rev.Lett. Ill (2013) 091302, [arXiv : 1304. 7791]. 

[82] V. C. Bust!, R. F. L. Holanda, and C. Clarkson, Supernovae as probes of cosmic parameters: 
estimating the bias from under-dense lines of sight, JCAP 1311 (2013) 020, [arXiv: 1309.6540]. 

[83] K. Kainulainen and V. Marra, A new stochastic approach to cumulative weak lensing, Phys. 
Rev. D80 (2009) 123020, [arXiv: 0909.0822]. 

[84] K. Kainulainen and V. Marra, Accurate Modeling of Weak Lensing with the sGL Method, Phys. 
Rev. D83 (2011) 023009, [arXiv: 1011.0732]. 

[85] K. Kainulainen and V. Marra, Weak lensing observables in the halo model, Phys. Rev. D 84 
(Sept., 2011) 063004, [arXiv: 1101.4143]. 

[86] V. Marra, M. Quartin, and L. Amendola, Accurate weak lensing of standard candles. 1. Flexible 
cosmological fits, Phys. Rev. D 88 (Sept., 2013) 063004, [arXiv: 1304.7689]. 

[87] M. Quartin, V. Marra, and L. Amendola, Accurate weak lensing of standard candles. II. 
Measuring erg with supernovae, Phys. Rev. D 89 (Jan., 2014) 023009, [arXiv: 1307.1155]. 

[88] L. Amendola, T. Castro, V. Marra, and M. Quartin, Constraining the growth of perturbations 
with lensing of supernovae, MNRAS 449 (May, 2015) 2845-2852, [arXiv: 1412.3703]. 

[89] L. Amendola, K. Kainulainen, V. Marra, and M. Quartin, Large-scale inhomogeneities may 
improve the cosmic concordance of supernovae, Phys. Rev. Lett. 105 (2010) 121302, 

[arXiv: 1002. 1232]. 

[90] K.-D. Nguyen Thu Lam and J. Kurchan, Stochastic perturbation of integrable systems: a 
window to weakly chaotic systems, ArXiv e-prints (May, 2013) [arXiv: 1305.4503]. 

[91] B. I. Halperin, Creen’s Functions for a Particle in a One-Dimensional Random Potential, 
Physical Review 139 (July, 1965) 104-117. 

[92] V. Perlick, Cravitational Lensing from a Spacetime Perspective, Living Reviews in Relativity 7 
(Sept., 2004) 9. 

[93] W. Feller, An introduction to probability theory and its applications. Wiley Series in Probability 
and Mathematical Statistics. New York etc.: John Wiley and Sons, 1971. 

[94] L. Arnold, Stochastic Differential Equations: Theory and Applications. Wiley, 1974. 

[95] I. Karatzas and S. E. Shreve, Brownian Motion and Stochastic Calculus. Springer-Verlag 
Berlin/Heidelberg/New York, 1991. 

[96] R. S. Liptser and A. Shiryaev, Statistics of random processes. 1: General theory. 

Springer-Verlag Berlin/Heidelberg/New York, 2001. 

[97] B. Oksendal, Stochastie Differential Equations: An Introduction with Applications. 

Springer-Verlag Berlin/Heidelberg/New York, 2003. 

[98] H. Risken, The Fokker-Planck equation. Methods of solution and applications. Springer Series in 
Synergetics, Berlin, New York: Springer, |cl989, 2nd ed., 1989. 

[99] H. McKean, Stochastic Integrals. Academic Press, 1969. 


- 46 - 


[100] P. Fleury, H. Dupuy, and J.-P. Uzan, Interpretation of the Hubble diagram in a 
nonhomogeneous universe, Phys.Rev. D87 (2013), no. 12 123526, [arXiv: 1302.5308]. 

[101] Z. Stuchlik, An Einstein-Strauss-de Sitter model of the universe, Bui. Astronom. Inst. 
Czechoslovakia 35 (1984) 205. 

[102] J.-P. Uzan, G. F. R. Ellis, and J. Larena, A two-mass expanding exact space-time solution, Gen. 
Rel. Grav. 43 (2011) 191-205, [arXiv: 1005. 1809]. 

[103] M. Mars, F. C. Mena, and R. Vera, Review on exact and perturbative deformations of the 
Einstein-Straus model: uniqueness and rigidity results. General Relativity and Gravitation 45 
(Nov., 2013) 2143-2173, [arXiv: 1307.4371]. 

[104] G. Lemaitre, L’Univers en expansion, Ann. Soc. Sci. Bruxelles, A 53 (1933) 51. 

[105] A. Krasinski, Inhomogeneous Cosmological Models. Cambridge University Press, 1997. 

[106] P. Szekeres, Quasispherical gravitational collapse, Phys. Rev. D 12 (1975) 2941. 

[107] G. Darmois, Les equations de la gravitation einsteinienne. Memorial des sciences 
mathematiques, 1927. 

[108] W. Israel, Singular hypersurfaces and thin shells in general relativity, Nuovo Cimento B Serie 
44 (July, 1966) 1-14. 

[109] W. Israel, Singular hypersurfaces and thin shells in general relativity, Nuovo Cimento B Serie 
48 (Apr., 1967) 463-463. 

[110] B. Panter, A. F. Heavens, and R. Jimenez, The mass function of the stellar component of 
galaxies in the Sloan Digital Sky Survey, MNRAS 355 (Dec., 2004) 764-768, 

[astro-ph/0406299]. 

[111] M. Bartelmann and P. Schneider, Weak gravitational lensing, Phys. Rept. 340 (2001) 291-472, 
[astro-ph/9912508]. 

[112] P. Schneider, J. Ehlers, and E. E. Falco, Gravitational Lenses. Springer-Verlag 
Berlin/Heidelberg/New York, 1992. 

[113] N. Deruelle and J.-P. Uzan, Theories de la relativite. Belin, Paris, 2014. 

[114] M. Sasaki, Cosmological gravitational lens equation: Its validity and limitation, Proq. Theor. 
Phys. 90 (1993) 753-781. 

[115] R. K. Sachs, Gravitational waves in general relativity. 6. The outgoing radiation condition, Proc. 
Roy. Soc. Land. A264 (1961) 309-338. 

[116] S. Seitz, P. Schneider, and J. Ehlers, Light propagation in arbitrary space-times and the 
gravitational lens approximation, Class. Quant. Grav. 11 (1994) 2345-2374, 

[astro-ph/9403056]. 

[117] J. Etherington, On the definition of distance in general relativity, Phil. Mag. 15 (1933) 761. 

[118] P. Fleury, Light propagation in inhomogeneous and anisotropic cosmologies. PhD thesis, 
Universite Pierre et Marie Curie, Paris 6, Nov., 2015. arXiv: 1511.03702. 

[119] P. Fleury, C. Pitrou, and J.-P. Uzan, Light propagation in a homogeneous and anisotropic 
universe, Phys. Rev. D91 (2015), no. 4 043511, [arXiv : 1410. 8473]. 


- 47 - 


